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i  Introduction  and  General  Goals 

This  final  technical  report  contains  a  summary  of  the  activities  supported  under  the  Air 
Force  AFOSR  PRET  Grant  F49620-96-1-0329  during  the  period  1  July  1996  through  31 
December  2001.  This  project  was  concerned  with  an  interdisciplinary  research  program  to 
develop  and  test  new  computational  methods  for  optimization  based  design  and  sensitivity 
analysis  of  aerospace  systems.  The  program  integrated  scientific  and  computational  tools 
developed  at  AeroSoft,  Beam  Technologies,  Boeing  Defense  and  Space  Group 
with  new  sensitivity  techniques  and  optimization  algorithms  developed  at  Virginia  Tech 
and  Cornell.  The  "focus  of  the  program  was  fundamental  research  in  sensitivity  and  adjoint 
based  methods  for  design,  control,  and  optimization  of  complex  aerospace  systems  governed 
by  partial  differential  equations.  In  addition,  the  program  was  structured  to  promote  t  e 
transition  of  this  basic  research  to  Air  Force  laboratories  and  to  industry.  The  principal 
investigator  was  Dr.  John  A.  Burns.  Hatcher  Professor  of  Mathematics  at  Virginia  Tech. 


2  Introduction  to  the  Center  and  General 
Research  Goals 

The  Air  Force  Center  for  Optimal  Design  and  Control  (CODAC)  was  established  in  May 
1993  under  an  Air  Force  AFOSR  URI  Grant.  CODAC  is  an  interdisciplinary  research 
center  with  core  academic  participants  at  Virginia  Tech.  CODAC  has  made  significant 
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progress  in  the  areas  of  optimal  design  and  control  of  systems  governed  by  partial  dif¬ 
ferential  equations.  In  addition,  CODAC  has  an  established  record  of  interactions  with 
industry  and  Air  Force  laboratories.  The  current  industrial  partners  include  AeroSoft  and 
Boeing  Defense  and  Space  Group.  This  group  is  composed  of  a  proven  team  of  scientists 
and  engineers  from  small  high-tech  firms  and  major  aerospace  companies.  In  addition  we 
have  interactions  with  North  Carolina  State  University  and  Cornell  University.  These  uni¬ 
versity  groups  furnishes  valuable  expertise  in  optimal  design,  control  and  optimization  of  a 
variety  of  application  areas.  The  research  group  at  Virginia  Tech  has  been  at  the  forefront 
of  the  development  of  sensitivity  methods  for  optimal  design,  with  applications  to  shape 
optimization  for  fluid  flow  management,  and  provides  expertise  in  computational  methods 
for  optimization  and  for  simulation  of  fluid  dynamics.  The  industrial  partners  themselves 
contribute  an  immense  expertise  in  computational  fluid  dynamics,  finite  element  modeling 
and  computational  mechanics. 

The  effort  supported  under  AFOSR  PRET  Grant  F49620-96- 1-0329  was  built  on  a  highly 
integrated  interdisciplinary  research  program  with  three  primary  goals: 

•  To  develop  and  investigate  new  mathematical  and  computational  methods  for  sensitiv¬ 
ity  analysis  with  applications  to  optimal  design  of  those  aerospace  systems  described 
by  nonlinear  partial  differential  equations.  These  applications  include  (but  are  not 
limited  to)  shape  optimization  and  design  for  flow  management,  materials  processing, 
manufacturing,  combustion  and  high  speed  flows.  By  working  jointly  with  industrial 
partners  and  applying  the  results  to  real-world  Air  Force  problems  we  established  a 
direct  and  rapid  track  for  transitioning  new  techniques  and  software  to  industry. 

•  The  development  of  a  mathematical  foundation  for  the  construction  of  new  models 
and  the  investigation  of  the  relationships  between  these  models,  simulation  meth- 
ods,  sensitivity  analysis  and  optimization  algorithms.  The  objective  is  to  provide  a 
theoretical  framework  for  the  rigorous  analysis  of  design  algorithms  that  combine  nu¬ 
merical  simulation  codes,  approximate  sensitivity  calculations  and  optimization  codes. 
In  particular,  it  is  important  to  determine  under  what  conditions  a  given  numerical 
simulation  scheme  can  be  combined  with  a  specific  optimization  method  to  produce 
a  convergent  design  algorithm. 

•  The  development  of  a  computational  environment  and  high  level  software  tools  that 
engineers  can  effectively  use  to  design  and  optimize  aerospace  systems.  The  objective 
is  to  initially  provide  a  common  software  substrate  to  all  university  and  industrial 
partners  in  the  center  and  then  develop  this  toolbox  into  a  computational  environment 
for  design.  The  existence  of  this  software  tool  will  enable  rapid  dissemination  of  the 
technology  developed  in  this  project  and  help  facilitate  interactions  between  university 
researchers  and  industrial  partners.  The  outcomes  envisioned  are  new  and  practical 
computational  tools  for  use  in  a  wide  range  of  aerospace  design  problems. 

3  CODAC  Organization  and  Facilities 

Dr.  John  A.  Bums  is  the  Director  of  CODAC  and  is  responsible  for  the  day-to-day  opera¬ 
tions  of  the  Center  and  for  implementing  and  coordinating  the  research,  labor atory/mdustry 
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interactions,  and  educational  programs.  Dr.  Eugene  M.  Cliff  is  the  Director  for  Engi¬ 
neering.  CODAC  is  located  within  the  Interdisciplinary  Center  for  Applied  Mathematics 
(ICAM)  at  Virginia  Tech.  Dr.  Terry  Herdman  is  the  Director  of  ICAM.  The  Executive 
Advisory  Committee  consists  of  the  Center  Director,  the  Director  for  Engineering  and  the 
Director  ICAM.  The  Executive  Committee  is  responsible  for  program  planning  and  for 
advising  the  Director  on  the  allocation  of  resources  and  on  ways  to  make  CODAC  more 
effective  as  an  Air  Force  resource.  The  research  conducted  at  CODAC  requires  computa¬ 
tional,  as  well  as  more  traditional  laboratory  facilities.  Although  much  of  the  initial  large 
scale  computing  was  done  on  supercomputers  at  Air  Force  facilities,  pre-processing  and 
post-processing  was  done  locally.  During  the  past  few  years  we  have  acquired  an  Origin 
2000  supercomputer  which  enables  us  to  develop  practical  software  tools  here  at  ICAhl. 


PRET  Center  Administration  and  Research  Team 

Research  under  the  PRET  Grant  was  conducted  at  six  institutions:  Virginia  Tech, 
Cornell,  AeroSoft,  Beam  Technologies,  Boeing  and  Lockheed  Martin.  Frequent  meetings 
of  the  team,  exchange  of  visits,  sharing  of  software,  exchange  of  graduate  students  and 
postdocs,  an  annual  industry- Air  Force  laboratory-university  workshop  and  communication 
of  results  were  be  coordinated  by  the  AFOSR-sponsored  PRET  Center  CODAC.  Professor 
John  Burns,  headed  an  advisory  board  which  consisted  of  Walters  (AeroSoft),  Roetman 
(Boeing)  and  Grossman  (MADD/Virginia  Tech).  Dr.  Grossman  was  the  Director  of  the 
Multidisciplinary  Analysis  and  Design  Center  at  Virginia  Tech,  has  broad  contact  with  the 
industrial  aerospace  design  community. 

The  PRET  Center  provides  the  structure  for  unique  educational  and  scientific  training 
of  students  and  post-doctoral  researchers  in  an  interdisciplinary  team  approach  to  scientific 
and  engineering  research.  CODAC  provided  unique  opportunities  for  theoretical,  compu¬ 
tational,  and  experimental  research.  Through  the  interactions  with  Air  Force  laboratories 
and  industrial  partners,  students  are  exposed  to  real  problems.  The  combined  theoretical, 
computational,  and  experimental  approach  provides  a  meaningful  interdisciplinary  research 
experience. 


ICAM  Computing  Facilities 

ICAM  houses  a  heterogeneous  Unix  system  with  file-sharing  under  a  Network  File 
System  (NFS)  and  has  excellent  computational  facilities.  This  includes: 

•  a  32-processor  SGI  Origin  2000  supercomputer  with  over  20  GB  of  RAM  for  large 
scale  computing, 

•  a  4-processor  SGI  Origin  200  server  providing  NFS  and  NIS  services, 

•  thirteen  Unix  workstations  [SGI  (9),  DEC  Alpha  (2),  Sun  (2)], 

•  a  collection  of  Pentium-based  PC’s  (5)  and  Power-PC  Macintosh  systems  (4)  with 
network  access  to  the  Unix  systems  and  the  Internet, 
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•  a  local  (OC3)  ATM  network  for  NFS  services  and  a  (100  BaseT)  switched  Ethernet 
network  for  user  services. 

The  research  developed  under  this  contract  has  wide  applicability  and  promises  con¬ 
siderable  payoff  in  aerospace  design  applications.  Although  we  worked  with  a  variety  of 
industrial  groups,  the  interaction  with  AeroSoft  Inc.  was  clear  the  most  fruitful.  This 
interaction  will  be  detailed  in  the  sections  below. 

lead  to  In  order  to  provide  focus  for  the  research  and  to  expedite  its  transition  to 
industrial  use  we  have  developed  research  partnerships  with  the  following  groups: 

Objectives 

The  research  program  had  three  primary  objectives:  (i)  To  develop  and  investigate  new 
mathematical  and  computational  methods  for  sensitivity  analysis  with  applications  to  op¬ 
timal  design  and  control  of  aerospace  systems;  (ii)  To  develop  a  mathematical  foundation 
for  the  construction  of  new  models  and  to  investigate  the  relationships  between  these  mod¬ 
els,  simulation  methods,  sensitivity  analysis  and  optimization  algorithms;  (iii)  To  develop  a 
computational  environment  and  high  level  software  tools  that  engineers  can  effectively  use 
to  design  and  optimize  aerospace  systems.  The  goal  was  to  provide  a  theoretical  framework 
for  the  rigorous  analysis  of  design  algorithms  that  combine  numerical  simulation  codes, 
approximate  sensitivity  calculations  and  optimization  codes. 

4  Accomplishments 

Under  the  support  of  this  grant  we  completed  the  development  a  suite  of  computational 
tools  for  sensitivity  analysis,  design,  control,  and  optimization  of  a  wide  variety  of  nonlinear 
partial  differential  equations.  We  focused  on  the  refinement  of  the  hybrid  continuous  sensi¬ 
tivity  methods  that  now  form  the  foundations  for  the  AeroSoft  product  SENSE  and  Beam 
Technologies’  CFETools.  We  made  considerable  progress  on  algorithms  for  coupling  single 
disciplinary  sensitivity  codes  for  application  to  multi-physics  problems.  This  research  has 
already  been  applied  to  COIL  Laser  design.  In  addition,  we  have  completed  an  analysis  of 
the  role  that  sensitivities  play  in  time  marching  numerical  schemes  and  developed  a  new 
fast  algorithm  for  solving  optimal  control  problems  governed  by  partial  differential  equa¬ 
tions  that  arise  in  fluid  flows.  This  work  was  motivated  by  interactions  with  our  Boeing 
partners  and  the  Air  Vehicles  Directorate  of  AFRL.  Although  this  grant  produced  several 
major  breakthroughs  in  sensitivity  computations  and  optimal  design,  we  will  focus  only  on 
three  major  projects.  These  projects  are  summarized  in  the  accomplishment  section  below. 
In  addition,  during  this  period  the  PRET  Center  has: 


•  Produced  more  than  115  scientific  papers,  articles  and  books, 

•  Made  more  than  150  presentations  at  conferences  and  colloquium, 

•  Directed  more  than  20  graduate  students  and  7  postdocs, 
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•  Worked  with  more  than  50  visitors,  representing  10  different  countries, 

•  Assisted  in  the  development  of  AeroSoft’s  SENSE  software  package, 

•  Worked  with  several  Air  Force  laboratories  to  transition  the  research  into  problems 
ranging  from  wind  tunnel  design  to  combustion  control. 


Major  Accomplishments 

We  give  a  brief  description  of  some  research  accomplishments  and  provide  an  indication  of 
the  significance  and  potential  applications,  details  are  given  in  the  next  section. 

Accomplishment  1:  Developed  a  Method  for  Estimating  Stability  Derivatives. 

•  New  Findings:  This  extends  work  begun  last  year  to  the  case  of  three-dimensional 
and  viscous  flows.  The  sensitivity  equation  provides  an  efficient  way  to  compute  static- 
stability  derivatives  which  describe,  in  linear  approximation,  the  way  the  aerodynamic 
forces  and  moments  depend  on  the  orientation  of  the  body  to  the  free-stream  flow. 
This  research  is  discussed  in  §  ??. 

•  Significance  and  Potential  Applications:  While  deducing  some  of  these  stability 
derivatives  using  finite-differences  of  neighboring  solutions  is  possible,  such  methods 
may  not  be  efficient.  For  example,  in  order  to  evaluate  stability  derivatives  associ¬ 
ated  with  asymmetric  flight,  nonlinear  flow-solutions  must  be  generated  on  an  entire 
configuration,  even  in  the  usual  case  wherein  the  vehicle  has  a  symmetry  plane.  The 
present  approach  allows  one  to  compute  derivatives,  such  as  the  weathercock  stability 
parameter  (Cn0  |/3=o),  based  on  a  nonlinear  flow  solution  computed  on  one-half  of  the 
configuration. 


Accomplishment  2:  Discovered  an  Extreme  Sensitivity  in  Nonlinear  Boundary 
Value  Problems  of  the  Type  that  Describe  Fluid  Flows. 

•  New  Findings:  For  convection  diffusion  problems  we  show  that,  because  of  finite 
precision  arithmetic  inherent  in  digital  computers,  convergent  numerical  algorithms 
can  produce  false  (purely  numerical)  solutions.  It  is  shown  that  these  false  solutions 
may  be  viewed  as  solutions  of  nearby  problems  with  very  high  sensitivity  to  boundary 
values.  In  addition,  because  of  this  high  sensitivity,  mesh  refinement  to  “converge  a 
solution”  only  exacerbates  the  matter. 

•  Significance  and  Potential  Applications:  The  above  results  show  that  “numerical 
based  proofs”  of  non-uniqueness  of  solutions  to  hydrodynamic  equations  (e.g.  Euler 
Equations)  need  further  study.  In  particular,  even  convergent  numerical  schemes  can 
produce  false  (non-unique)  numerical  solutions  that  are  not  solutions  of  the  underlying 
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boundary  value  problem.  This  is  extremely  important  when  such  codes  are  used  for 
optimal  design  and  control.  This  analysis  also  can  be  applied  to  the  development  of 
new  accelerated  marching  algorithms.  Finally,  the  sensitivity  and  stability  analysis 
developed  in  this  effort  can  be  applied  to  more  complex  fluid  flow  phenomena.  In 
particular,  it  supports  the  work  of  Bamieh,  Butler,  Dahleh,  Farrell  and  Trefethen  that 
suggest  transition  to  turbulence  in  certain  fluid  flows  are  a  result  of  flow  sensitivity. 


5  Technical  Details 

Here  we  present  a  few  technical  details  concerning  the  accomplishments  given  above.  These 
results  may  be  found  in  the  papers  listed  at  the  end  of  the  report. 

5.1  Direct  Calculation  of  Aerodynamic  Force  Derivatives:  A  Sensitivity- 
Equation  Approach 

In  this  section  we  discuss  the  sensitivity-equation  approach  to  computing  stability  deriva¬ 
tives  using  a  single  non-linear  solution  to  the  underlying  fluid  equations.  The  sensitivity 
equations  are  presented  in  integral  form  with  the  necessary  boundary  conditions.  The 
lift-curve  slope  is  computed  at  several  angles  of  attack  for  a  laminar  airfoil.  Stability  char¬ 
acteristics  are  analyzed  for  the  YB-49  flying  wing. 

5.1.1  Introduction 

Computational  Fluid  Dynamics  (CFD)  plays  an  increasingly  important  role  in  the  analysis 
and  design  of  aerospace  vehicles.  From  a  flight  mechanics  view  a  primary  purpose  of  CFD 
analyses  is  timely  and  accurate  prediction  of  the  aero-propulsion  forces  and  moments  applied 
to  the  vehicle.  For  control-system  design  we  need  to  have  some  notion  of  how  these  forces 
and  moments  vary  with  vehicle  motion  -  that  is,  one  needs  to  estimate  certain  stability 

derivatives.  . 

While  deducing  some  of  these  stability  derivatives  using  finite-differences  of  neighboring 
solutions  is  possible,  such  methods  may  not  be  efficient.  For  example,  in  order  to  evaluate 
stability  derivatives  associated  with  asymmetric  flight,  nonlinear  flow-solutions  must  be 
generated  on  an  entire  configuration,  even  in  the  usual  case  wherein  the  vehicle  has  a 
symmetry  plane.  The  present  approach  allows  one  to  compute  derivatives,  such  as  the 
weathercock  stability  parameter  (Cn0  \p=0),  based  on  a  nonlinear  flow  solution  computed 
on  one-half  of  the  configuration.  In  addition,  certain  derivatives  (e.g.  the  damping-m-roll 
parameter  Ctp),  require  a  flow  solution  that  permits  explicit  motions  of  the  body  (in  this 
case  roll-rate  p).  Such  modeling  may  not  be  readily  available  in  standard  CFD  codes. 

In  recent  years  the  desire  to  use  optimization-based  design  algorithms  has  spurred  the 
need  for  design  sensitivities  and  for  efficient  ways  to  calculate  them.  At  this  point  there  are 
three  approaches  to  the  calculation  of  design  sensitivities. 

1.  finite  difference  of  neighboring  solutions 
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2.  ‘differentiate’  the  numerical  code  (  i.e.,  ADIFOR) 

3.  ‘differentiate’  the  boundary-value  problem 

While  there  are  subtleties,  in  the  second  approach  formally  one  imagines  applying  the  chain- 
rule  to  each  line  of  code  to  produce  another  code  which  will  evaluate  the  derivatives.  Here 
we  will  focus  on  the  third  approach,  often  called  the  Sensitivity  Equation  Method  (SEM). 
A  high-level  view  of  the  abstract  SEM  approach  is  that  it  produces  a  linear  boundary- 
value  problem;  the  solution  of  this  problem  is  a  flow  sensitivity  and  describes,  in  linear 
approximation,  the  way  the  flow  solution  (dependent  variables)  depends  on  a  (scalar)  design 
parameter.  In  applications  one  finds  that  the  solution  of  the  underlying  nonlinear  boundary- 
value  problem  often  relies  on  a  linearization  so  that  many  of  the  required  ingredients  for 
the  SEM  already  exist  within  the  basic  CFD  code. 

A  related  approach  applies  the  SEM  to  the  discretized  equations.  This  is  also  related 
to  the  ADIFOR  approach,  in  that  both  are  applied  to  a  discretized  problem.  In  general  a 
discrete  approximation  of  the  linear  SEM  is  not  the  same  as  applying  sensitivity  ideas  to 
the  discrete  equations.  For  one  thing,  the  latter  approach  requires  grid  sensitivities  which 
need  not  arise  in  the  abstract  method.  There  are  connections  between  these  approaches, 
including  a  notion  of  asymptotic  consistency  develped  under  the  present  research. 

In  this  section  we  will  develop  and  demonstrate  the  abstract  SEM  method  for  computing 
stability  derivatives.  Thus,  we  must  formulate  the  flow  equations  so  that  the  desired  flight 
parameters,  such  as  angle-of-attack  and  side-slip  angle,  are  available  explicitly.  Thus,  the 
following  discussion  is  devoted  to  the  parameterization  of  the  body  geometry  and  to  a  careful 
description  of  some  underlying  coordinate  systems.  Following  this  we  will  briefly  describe 
the  flow-equations  and  the  associated  boundary  conditions.  The  main  contribution  is  the 
derivation  of  sensitivity  equations  for  this  class  of  problems.  We  demonstrate  the  approach 
with  several  examples. 

5.1.2  Body  Geometry 

In  CFD  applications  the  problem  geometry  is  commonly  described  in  terms  of  a  grid.  Ab¬ 
stractly,  one  imagines  a  map  from  a  computational  domain  to  the  physical  domain  and  the 
grid  provides  a  discrete  version  of  this  map.  The  body  geometry  is  implicitly  described 
by  imposing  appropriate  boundary  conditions  on  specific  planes  or  parts  of  planes  in  the 
computational  domain. 

For  our  purposes  we  must  describe  how  the  body  moves  and/ or  deforms  under  specific 
changes  in  the  geometry.  Since  we  are  concerned  with  relative  motions  between  the  body 
and  the  surrounding  fluid,  there  are  inevitably  several  reasonable  choices  for  underlying 
coordinates.  Our  approach  is  arguably  natural  from  a  flight  mechanics  point-of-view. 

We  assume  the  atmosphere  is  homogeneous  and  at  rest  ’and  use  the  symbol  Vc  to  denote 
the  lineal  velocity  of  a  specified  point  on  the  vehicle.  At  this  specified  point  we  imagine  two 
coordinate  systems  with  a  common  origin 

•  a  local  horizontal  system  OxhyhZh,  and 

•  a  body-fixed  system  Ox^y^. 
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The  local  horizontal  system  has  the  OxhVh  horizontal  with^  vertical  and  down.  Since 
the  atmosphere  is  homogeneous  and  at  rest  we  can  take  Vc  =  [1,  0,  0]£,  so  that  the 
basic  translational  motion  is  in  the  Xh  direction.  Because  of  this  choice,  the  body-fixed 
system  Oxbybzb  is  oriented  to  the  horizontal  system  by  the  aerodynamic  angles:  a,  the 
angle-of-attack;  and  /?,  the  sideslip  angle. 

The  jig-shape  of  the  vehicle  is  defined  by  a  map  n0  :  E  C  1R2  IR.3.  Here  a  €  £  are 

points  in  the  computational  plane  that  define  the  body  boundary  or  a  part  of  the  body.  In 
applications,  one  does  not  deal  with  the  map  7r0  but  rather  with  the  image  points  produced 
:  that  is,  the  points  of  the  physical  domain  that  correspond  to  appropriate  boundary  points 
in  the  computational  domain.  In  applying  the  abstract  SEM,  we  will  need  to  compute 
various  partial  derivatives  so  that  we  must  describe  explicitly  the  underlying  maps. 

We  are  concerned  with  changes  in  the  flow  field  (and  so,  changes  in  the  applied  loads) 
that  are  induced  by  deformations  and  motions  of  the  body.  For  body  deformations  we  admit 
a  deformation  field  tti  so  that  points  on  the  surface  of  the  deformed  body  are  given  by 

rs(cr,<£  i)=7f0(a)  +  (5-l) 

where  the  scalar  <j>\  is  the  magnitude  of  the  deformation.  It  is  possible  to  admit  a  com¬ 
bination  of  deformations  7Ti,7T2,  . . .,  but  since  we  deal  with  these  one  at  a  time  the  single 
deformation  case  is  adequate.  In  design  applications  the  field  tvi  would  describe  a  per¬ 
turbation  of  interest,  for  example,  geometric  twist  distribution  in  a  wing.  In  aeroelastic 
applications,  the  field  would  describe  some  elastic  mode  -  the  SEM  then  gives  a  direct 
calculation  of  the  associated  aeroelastic  stiffness  parameter. 

5.1.3  Integral  Equations  for  Fluid  Dynamics 

The  scope  of  this  paper  encompasses  both  inviscid  and  viscous  flow  fields.  Solutions  to  the 
fluid-dynamic  equations  provide  coefficients  for  the  sensitivity  equations  which  are  discussed 
below.  The  three-dimensional  flow  of  a  calorically-perfect,  viscous  fluid  is  governed  by 
a  system  of  non-linear,  hyperbolic  partial  differential  equations  which  can  be  written  in 
integral  form  as 

^  III  ^  ^  (F(Q)  •  n)  dA  =  j)  (F„(Q)  •  n)  dA.  (5.2) 

The  conservative  variables  are  Q  =  Q  (x,y,z,t)  =  [p,pu,pv,pw,peo]T  and  represent  the 
density,  momentum  and  total  energy  per  unit  volume  of  the  fluid.  The  surface  integrals 
represent  the  inviscid  and  viscous  fluxes  (F  and  F„).  Transport  properties  are  computed 
using  Sutherland  curve  fits.  The  integral  formulation  is  the  fundamental  starting  place  for 
finite- volume  discretizations. 

This  system  of  equations  may  be  solved  subject  to  boundary  conditions  which  for  an 
inviscid  flow  involve  scalar  relationships  such  as 

Vn  =  0  (5-3) 

V/n  =  0,  (5-4) 

where  the  velocity  vector  is  depicted  as  V  =  (u,  v,  w)T .  For  a  viscous  fluid,  the  no-slip 
condition  is  simply  V  =  0. 
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While  the  model  Eqn.  (5.2)  includes  explicit  time-dependence,  we  are  interested  only 
in  the  steady-state  solution,  which  is  assumed  to  exist.  For  this  purpose  we  introduce  the 
notation 


Q (x,y.z)  =  lim  Q (x,y,z,t). 

t— +OO 


5.1.4  Sensitivity  Equations 

The  idea  of  a  flow  sensitivity  is  that  of  a  linear  approximation  -  loosely  a  partial  derivative 
of  the  flow  with  respect  to  a  parameter  of  interest.  Here  we  use  the  symbol  rj  to  denote 
a  generic  parameter  and  to  emphasize  that  the  flow  solution  depends  both  on  position  in 
space  and  on  the  parameter  we  write 

Q  (x,y,z;  rj). 


The  sensitivity  we  seek  is  then  formally  given  by 


S  = 


5Q 

dr} 


(5.5) 


Our  objective  is  to  derive  a  linear  boundary- value  problem  for  S(-;  77).  We  do  this  formally, 
by  differentiating  the  elements  of  the  nonlinear  boundary- value  problem  for  the  flow  Q(-;  rj). 
Here  we  are  specifically  interested  in  two  parameters: 


•  a  -  the  angle-of-attack,  and 


•  /3  -  the  side-slip  angle. 


5.1.5  Finite- Volume  Sensitivity 


Our  form  of  the  governing  equations  in  Eqn.  (5.2)  is  valid  in  an  inertial  reference  frame.  In 
particular,  we  choose  a  frame  such  that  the  Cartesian-coordinate  axes  agree  with  with  the 
instantaneous  body  axes  (jig  shape).  The  parameters  of  interest  (a,  0)  then  do  not  explicitly 
appear  in  any  of  the  flux  functions.  In  this  case  we  proceed  by  formally  differentiating 
Eqn.  (5.2)  with  respect  to  77  and  then  interchange  orders  of  differentiation  to  find 

jt  JJJ  S  dV  +  £  (f'(Q,  S)  •  n)  dA  =  jf  (f^Q,  S)  •  n)  dA,  (5.6) 


where  Q  represents  the  conservative  variables  available  from  the  steady-flow  solution  and 
S  represents  the  unknown  flow'  sensitivities.  Note  that  the  governing  equations  for  S(-;  77) 
is  linear  and  that  the  flow  solution  ( Q(-,r ?)  )  enters  through  spatially  varying  coefficients. 
We  have  previously  demonstrated  that  in  the  Euler  flow  case,  Eqn.  (5.6)  is  equivalent  to 
computing  the  77-derivative  of  the  differential  form  of  the  conservation  law. 

The  inviscid  flux  vector  in  the  sensitivity  equations  is  written  as 


F'  n 


dF  dQ  .  0FO  . 
<9Q  dr)  <9Q 


< 


v. 


p'(V-n)  +  p(V'.n) 

p'u(V  •  n)  +  /9t/(V  •  n)  +  pu(V'  •  n)  +  hxp' 
p'v{V  ■  n)  +  pv'(V  ■  n)  +  pv(V^  ■  n)  +  nyp' 
p'w(V  ■  n)  +  pw'iy  ■  n)  +  pw{W'  ■  n)  +  hzp' 
p'ho(V  •  n)  +  ph'0(V  ■  n)  +  ph0(V'  •  n) 


}  ,  (5-7) 
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where  the  density  sensitivity  is  p',  the  velocity  sensitivity  is  V'  =  (u1,  v',wr)T,  and  the 
pressure  sensitivity  is  p'.  The  stagnation  enthalpy  per  unit  mass  is  defined  as  ho  =  h  +  (V  • 
V)/2  and  the  enthalpy  is  h  =  e  +  p/p.  The  internal  energy  for  a  calorically  perfect  gas  is 
e  =  5/2  RT.  The  numerical  flux  is  computed  using  a  characteristic-based,  upwind  scheme. 

The  viscous  flux  vector  may  be  written  as 


F'vh={ 


0 


rxxnx  “h  ^xy^y  ^xz^z 

ryx™x  Tyy™y  Tyz™z 

r'zxnx  +  T'zyhy  +  t'z 


n. 


(5.8) 


Assuming  a  linear  stress-strain  relationship,  the  viscous  shear  can  be  written  as  Tjj  — 
2 p(Sij  -  1/3  Uktk$ij)  where  S{j  represents  the  strain  rate.  The  sensitivity  of  the  shear- 
stress  tensor  is  derived  by  differentiation 


t'j  =  2p' 


?y-§(V*V)«y 


+  2/i 


(5.9) 


where  the  sensitivity  strain  rate  is 


The  sensitivity  of  the  conduction  term  is  ( q )'  =  —k'VT  —  kVT 


5.1.6  Boundary  Condition  Sensitivity 

Since  our  parameters  do  not  explicitly  enter  the  flux  terms  in  Eqn.  (5.2),  the  resulting 
sensitivity  equation  Eqn.  (5.6)  is  homogeneous.  All  the  action  then  is  in  the  boundary 
conditions  and  in  their  explicit  dependence  on  the  problem  parameters. 


In-Flow  Condition 

Our  model  is  expressed  in  a  coordinate  system  translating  with  the  vehicle  so  that  the 
relative  fluid  velocity  is  given  by  the  usual  Gallilean  transformation 

V  =  Vabs  —  Vc. 


In  the  present  case  all  of  the  fluid  at  the  in-flow  boundary  is  at  rest,  so  that  we  have  the 


inflow  conditions: 


v(*0  =  -l|vc|| 


cos  / 3  cos  a 
sin  (3 

cos/?  sin  a  < 


The  in-flow  boundary  condition  for  the  a-sensitivity  is  accordingly 


3q 

da 


-  liVcII  { 


0 

cos  (3  sin  a 

0 

—  cos  /?  cos  a 


(5.10) 
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while  for  the  /3-sensitivity  we  have: 


0 


<9q 

00 


=  l|Vc||  l 


sin  0  cos  a 
—  cos/3 
sin  /3  sin  a 
0 


>• 


(5.11) 


No-Penetration  Condition 

Another  of  the  important  boundary  conditions  in  Euler  flows  is  the  so-called  no  pene¬ 
tration  condition  (5.3),  now  written  as  V  •  ft  =  0  on  the  solid  boundary.  Using  the  notation 
introduced  above  we  can  write  this  more  precisely  as 


V (7r(<7,  77) ;  77)  *  n(cr,  77)  —  0  Vcr  €  £*>. 


(5.12) 


Since  (5.12)  must  hold  for  all  t]  in  some  open  region,  we  can  formally  differentiate  with 
respect  to  77.  Note  that  the  flow  velocity  function  is  defined  over  the  flow  region,  that  is 
V  :  O  x  Clv  R3,  where  O  C  R3  is  the  flow  region.  In  computing  the  77-derivative  of  (5.12) 
one  must  be  careful  to  account  for  the  chain-rule  terms.  This  leads  to 


(  dv  dix  av 
(  dx  dr]  +  dr] 


•n  + V- 


=  0. 


The  velocity  sensitivity  we  seek  is 


captured  in  the  ( 


term  so  that  we  are  led  to  define 


This  is  the  sensitivity  variable  we  are  seeking  and  the  boundary  condition  on  S  that  arises 
from  the  flow-tangency  condition  is 


S 


•  n  = 


avavr  .  -  (dh\ 

dx  dr]  \dr] ) 


(5.13) 


The  right-side  of  the  boundary  condition  Eqn.  (5.13)  can  be  interpreted  as  a  transpiration 
condition  for  the  sensitivity.  As  noted  above,  the  state  of  the  fluid  in  terms  of  primitive  flow 
variables  at  a  given  point  in  space  is  given  by  five  quantities  q  =  [p  u  v  w  p]T  In  equation 
(5.13)  we  are  describing  the  sensitivity  of  the  velocity  components.  Thus,  if  we  wish  to  use  S 
to  denote  the  five-component  state-sensitivity  there  is  a  3  x  5  matrix  projection  introduced. 
The  left  side  of  (5.13)  can  be  written  as: 

MS  ■  n  —  . . .  or  equivalently  S  •  (MTn)  =  . . . 

The  particular  parameters  of  interest  in  the  present  case  are  the  angle-of-attack  a  and  the 
side-slip  angle  0. 

Additional  Solid  Boundary  Conditions 
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Beyond  the  in-flow  and  flow-tangency  conditions  it  is  common  to  introduce  additional 
boundary-conditions  at  solid  interfaces  in  Euler  flows.  In  the  present  case  we  shall  use  three 
additional  (scalar)  conditions.  Each  of  these  is  written  in  a  standard  form  as 

Vi  [f(q(x; »?))]  •  n(x, rj)  =  0  VxeT  (5.14) 

where  T  denotes  the  boundary  points.  In  our  application  we  propose  to  use  the  following 
three  scalar-valued  functions  for  /. 

1.  fi(q)  =  qi  =  P  density 

2.  f2(q)  =  Qs=P  pressure 

3.  /3O  =  (<72 2  +  <73 2  +  <?42)/2  kinetic  energy 

Note  that  the  scalar-functions  (/)  are  not  explicitly  functions  of  the  design-parameter  (77). 
We  exploit  this  to  expand  (5.14)  as 

Vx9(tt(<7;  77)577)  -n(a,  77)  =  0  V  a  €  Eb. 

Recall  that  points  on  the  boundary  are  images  of  the  map  7r.  The  77  derivative  of  this  leads 

VXS  •  h  =  —  (§0  [(V«*  q)nv  •  n  +  V*9  •  A,] .  (5.15) 

The  first  of  the  forcing-terms  on  the  right  generally  involves  the  Hessian  of  the  primitive 
variables.  However,  for  the  present  application  with  77  ==  a  or  /3  we  have  that  both  7r,j  and  n-q 
vanish,  so  that  in  our  applications  the  b.c.  (5.15)  are  homogeneous. 

5.1.7  Applications 

As  a  first  numerical  example  of  the  SEM  approach  to  estimating  stability  derivatives,  we 
consider  the  calculation  of  the  lift-curve  slope  for  a  2-D  viscous,  compressible  flow  around  an 
airfoil.  A  potential  benefit  for  design  engineers  is  the  ability  to  determine  the  lift-curve  slope 
using  a  single  non-linear  flow  solution.  Specifically,  we  model  the  flow  around  a  NLF(l)- 
0416  airfoil  developed  at  NASA  Langley  by  Somers.  The  flow  conditions  (Mx>  =  0.5, 
Re  =  2  x  106)  are  extracted  from  an  AGARD  set  of  test  cases  for  validating  CFD  codes. 
The  pressure  contours  and  streamlines  at  a  =  0°  are  shown  in  Fig.  5.1  on  a  185  x  97  x  2 
“C”  mesh. 

To  proceed  in  the  analysis,  we  wish  determine  the  sensitivity  of  the  airfoil’s  lift  to  the 
angle  of  attack.  The  far-field  boundary  condition  for  the  S£M  is  determined  by  writing  the 
free-stream  velocity  as  a  function  of  the  magnitude  and  angle  of  attack.  In  the  absence  of 
side-slip,  we  apply  Eqn.  (5.10)  as  our  free-stream  sensitivity  vector  with  0  =  0°. 

Using  this  free-stream  boundary  condition,  the  solution  to  the  sensitivity  equations  at 
a  =  0°  is  shown  in  the  following  figure.  Sensitivity  pressure  contours  and  streamlines  of 
velocity-vector  sensitivity  are  depicted.  From  the  pressure  sensitivity  near  the  nose,  we  see 
that  an  increase  in  the  angle  of  attack  will  cause  an  increase  in  the  lower-surface  pressure 
and  a  decrease  in  the  upper-surface  pressure.  The  streamlines  show  that  the  magnitude  of 
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Figure  5.1:  Pressure  contours  and  velocity  vectors  for  the  NLF(1)-0416  airfoil  at  zero  angle 
of  attack. 
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the  velocity  vector  will  decrease  beneath  the  airfoil  surface  and  increase  along  the  upper 
surface.  Careful  investigation  at  the  trailing  edge  shows  that  the  momentum  of  the  boundary 
layer  is  predicted  to  decrease.  All  these  indications  are  consistent  with  an  increase  in  the 
angle  of  attack  for  attached  flow.  The  pressure  coefficient  and  sensitivity  are  given  in 
FigFig:NLF4i3. 


Figure  5.2:  Pressure  sensitivity  and  streamlines 


NLF(1)-0416  Pressure  Distribution 
M^O.5,  Re=2x10#,  a=0  deg. 


Figure  5.3:  Pressure  coefficient  and  sensitivity 

To  determine  the  accuracy  of  the  sensitivity  method  (more  fundamentally,  the  linearity 
of  the  local  flow  condition) ,  we  extrapolate  the  baseline  flow  to  a  =  4°  using  a  Taylor-series 
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expansion 


Q=QW+^ 


(5.16) 


W 

Aa. 


The  pressure  and  streamlines  for  the  near-by  approximation  is  given  in  Fig.  5.4.  The  pre¬ 
dominant  flow  features  are  a  rearward  movement  of  the  stagnation  point  at  the  nose  and 
decreased  pressure  along  the  upper  surface.  Notice  that  the  streamlines  do  not  indicate  sep¬ 
arated  flow.  The  accuracy  of  the  pressure  coefficient  on  the  airfoil  surface  can  be  evaluated 
using  Fig.  ??. 


Figure  5.4:  Projected  solution  at  a  =  4° 

Performing  a  sensitivity  analysis  at  a  =  4°  produces  the  contour  solution  shown  in 
Fig.  5.6.  The  bubbles  in  the  upper-surface  pressure  sensitivity  indicate  the  beginning 
and  ending  of  a  laminar  separation  bubble  followed  by  aerodynamic  stall.  The  sensitiv¬ 
ity  streamlines  depict  similar  behavior  to  the  zero-degree  case  beneath  the  airfoil.  However, 
the  prediction  for  the  flow  above  the  airfoil  is  for  further  flow  separation.  This  is  seen  by 
the  upward  movement  of  the  sensitivity  streamlines.  Assuming  laminar  flow  over  the  entire 
airfoil  surface,  the  computational  solution  begins  to  separate  at  a  =  5°,  consistent  with  the 
sensitivity  calculation. 

The  surface  pressure  coefficient  and  sensitivity  at  a  =  4°  is  shown  in  Fig.  5.7.  The 
sensitivity  profile  is  much  different  than  the  smooth  variation  from  nose  to  tail  shown  at 
a  =  0°.  The  sensitivity  further  indicates  the  likelihood  of  separated  flow  at  higher  angles  of 
attack.  The  Taylor  approximation  at  a  =  8°  shows  a  pressure  plateau  on  the  upper  surface 
which  indicates  the  presence  of  a  laminar  separation  bubble. 

Our  objective  of  determining  the  lift-curve  slope  from  one  flow  calculation  and  one 
sensitivity  calculation  can  be  achieved  by  differentiating  the  lift  coefficient  with  respect  to 
the  angle  of  attack.  Writing  the  lift  coefficient  in  terms  of  force  coefficients  in  the  x  and  y 
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NLF(1)-0416  Pressure  Distribution 

M^O.5,  Re=2x10*.  a=4  deg. 


Figure  5.5:  Surface  pressure  coefficient 


direction,  we  have 

Ci  =  Cy  cos  a  —  Cx  sin  a, 
so  that  the  lift-curve  slope  is  computed  as 

=  f  — 1M.  cos  a  —  dfi  sin  —  (Cy  sin  a  +  Cx  cos  a)  (5-17) 

da  \  da  da  J 

The  sensitivity  of  the  two  force  coefficients  is  obtained  by  integrating  the  inviscid  and  viscous 
fluxes  on  the  airfoil  surface.  The  lift-curve  slope  is  shown  in  Fig.  5.8  losing  experiment, 
computation  and  sensitivity  calculations.  The  computational  lift-curve  slope  shows  more 
lift  at  the  same  angle  of  attack  when  compared  to  experiment.  The  CFD  predicts  a  higher 
pressure  on  the  lower  surface,  thus  producing  more  lift.  The  sensitivity  lift-curve  slope  is 
given  at  a  =  0°  and  4°. 

5.1.8  Flying-Wing  Stability  Analysis 

As  early  as  1923,  Jack  Northrop  began  advocating  the  flying-wing  concept  as  the  next 
revolution  in  aircraft  design.  In  the  late  1930’s  the  Northrop  Corporation  undertook  the 
development  of  a  series  of  all-wing  concepts;  and  in  late  1941  received  an  order  for  two 
propeller-driven  flying  wings  in  support  of  the  long-range  bombing  requirements  of  the 
U.S.  Army  Air  Corps.  These  XB-35’s  were  plagued  with  both  maintenance  and  accident 
difficulties.  With  the  growth  of  the  jet-age,  the  XB-35  airframe  was  updated  with  eight 
Allison  J35-A-5  turbo-jets  under  a  contract  issued  in  June,  1945.  The  new  aircraft  was 
designated  the  YB-49  and  a  schematic  of  the  flying-wing  bomber  is  shown  in  Fig.  5.9. 
Airframe  parameters  for  the  plane  are  given  in  Table  1. 

The  perceived  advantages  of  the  flying  wing  were  increased  aerodynamic  efficiency  from 
reduced  parasite  drag,  simpler  construction  methods  from  fewer  structural  complications, 
and  increased  stealth  from  a  smaller  visual  target.  Unfortunately,  the  technology  of  the  day 
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0.5 


Figure  5.6:  Sensitivity  of  the  pressure  and  velocity  vectors  to  the  angle  of  attack  at  a  —  4°. 
Assuming  totally  laminar  flow,  the  airfoil  begins  to  stall  at  a  =  5°. 

could  not  solve  problems  associated  with  the  aircraft’s  stability  characteristics  -  namely, 
difficulty  holding  a  steady  course  and  altitude,  as  well  as  poor  damping  in  yaw.  The  fly-by¬ 
wire  concept  alleviates  some  of  these  problems  in  the  B-2  stealth  bomber. 

In  this  section,  we  will  use  the  sensitivity-equation  method  to  estimate  some  of  the 
stability  characteristics  of  the  YB-49  flying  wing.  To  simplify  grid  generation,  we  have 
neglected  the  bubble  canopy,  and  four  vert ical-fins/wing- fences  that  are  a  part  of  the  actual 
YB-49.  Our  calculations  sample  the  angle-of-attack  range  (a  =  0°,4o,8°)  and  we  compute 
CxoJ  CMa,  Cn0  and  Q0  derivatives.  The  free-stream  conditions  correspond  approximately 
to  cruise  speed  (M0 0  =  0.5)  at  an  altitude  of  35,000  feet. 

Pitch  Stability 

Stability  in  pitch  requires,  of  course,  that  in  the  neighborhood  of  an  equilibrium  point 
the  slope  of  the  pitch  moment  with  angle-of-attack  be  negative.  This  leads  to  the  stipulation 
that  the  center-of-gravity  be  suitably  forward  of  the  mean  aerodynamic  center  -  defined, 
loosely,  as  the  point  about  which  the  moment-slope  is  zero. 

For  an  inviscid-flow  CFD  model,  the  pitching  moment  is  computed  as  an  appropriate 
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NLF(1)-0416  Pressure  Distribution 
M^O.5,  Re=2x106,  Projection  Near  Stall 


Figure  5.7:  Surface  pressure  coefficient  at  a  -  4°  and  8°  along  with  the  pressure  sensitivity. 


surface  integral  of  the  normal  pressure.  We  write  this  as  an  iterated  integral 


M(a) 


7*2  (  ~2 

p(x,y,z\a )  d( — - — ) 


dy. 


(5.18) 


The  inner  term  is  identified  as  a  familiar  section  contour  integral  and  represents  the  con¬ 
tribution  to  the  pitching  moment  at  spanwise  location  y.  One  may  plot  this  pitch-loading 
against  the  spanwise  location  y  as  shown  in  Fig.  5.10. 

In  the  present  treatment  the  longitudinal  stability  parameter  Ma  is  evaluated  in  the 
same  way  with  the  pressure  replaced  by  the  pressure  sensitivity  pa.  One  may  describe  the 
inner  integral  as  the  stability  loading  function;  it  describes  the  contribution  to  pitch  stability 
from  various  locations  along  the  wing-span.  Fig.  5.11  displays  the  stability  loading  for  the 
YB-49  planform  at  a  =  4°  and  a  =  8°.  While  the  result  is  not  dramatic,  one  can  clearly 
see  the  decrement  in  the  stability  contribution  from  the  outboard  wing  panels.  This  feature 
was  described  qualitatively  by  J.  Northrop  in  his  1947  Wilbur  Wright  Memorial  Lecture. 

Contours  of  the  pressure  sensitivity  to  the  angle  of  attack  are  shown  in  Fig.  5.12  and 
Fig.  5.13.  The  domain  is  composed  of  ten  zones  composed  of  approximately  620,000  grid 
points.  Both  the  Euler  equations  and  the  sensitivity  equations  are  solved  using  a  first- 
order  spatially  accurate  flux-vector  splitting  scheme.  The  residual  for  the  Euler  equations 
are  converged  four  orders  of  magnitude  in  approximately  2250  iterations.  The  sensitivity 
equations  exhibit  some  stiffness  for  these  calculations  that  has  not  been  experienced  in 
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Lift  Curve  for  NLF(1)-0416 


Angle  of  Attack,  a  (deg) 


Figure  5.8:  Lift-curve  slopes  obtained  from  experiment,  CFD  flow  solver  and  the  sensitivity- 
equation  method. 


Figure  5.9:  Three- view  schematic  of  the  YB-49  flying  wing. 


past  SEM  calculations.  The  theoretical  reasons  for  this  have  not  been  investigated.  In 
general,  the  sensitivity  equations,  being  linear,  can  be  solved  using  an  infinite  time  step.  A 
smaller  time  step  is  used  for  the  77  =  or,  /3  cases  for  two  reasons:  the  fore-mentioned  stiffness 
and  also  because  the  problem  is  three-dimensional  which  requires  a  compromise  between 
convergence  and  memory.  A  hybrid  Euler  implicit/relaxation  scheme  is  used  to  sweep 
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Wing  span  (6) 
Wing  area  ( S ) 
Aspect  ratio  (A) 
Taper  ratio  (A) 

52.4  m 

371.6  m2 

7.4 

0.25 

Chord 

11.4  m  (Root) 

2.8  m  (Tip) 

Airfoil  section 

NACA  653  -  018  (Root) 
NACA  653  -  019  (Tip) 

Incidence 

0°  (Root) 

-4°  (Tip) 

0°53' 

Dihedral 

Sweep 

26°57'48”  (LE) 

10o15'22”  (TE) 

M.A.C. 

8  m  (7.33  m  aft  of  nose) 

C.G. 

7.13  m  aft  of  the  nose  (est.) 

Ceiling 

45,700  feet 

Cruising  speed 

429  mph  at  35,000  feet 

Table  1:  Airframe  data  for  the  YB-49  configuration. 


through  the  domain.  Updates  to  the  sensitivity  vector  are  computed  on  two-dimensional 
“planes”  (in  computational  space)  using  Jacobi  inner  iterations  with  non-linear  residual 
updates  computed  in  the  third  direction.  Therefore,  a  complete  linearization  of  the  SEM 
residual  requires  storage  of  the  entire  linearization  matrix  for  all  computational  planes  and 
all  zones  which  is  impractical  with  limited  memory.  The  compromise  made  for  saving 
memory  is  diminished  convergence  rate.  The  sensitivity  calculations  require  an  average  of 
60  iterations  to  converge  -  about  2.7%  of  the  time  to  compute  the  flow  solution.  The  SEM 
equations  requires  more  memory  than  the  Euler  solutions  because  both  the  flow  variables 
(Q)  and  sensitivity  variables  (Q')  are  retained  in  memory  during  the  calculation. 

The  performance  of  the  wing’s  lift  and  moment  are  summarized  in  Table  2  and  Table  3, 
respectively.  A  comparison  is  made  with  a  standard  semi-empirical  formula  for  the  lift- 
curve  slope  which  corrects  the  two-dimensional  airfoil  Cia  for  taper  ratio,  sweep  angle, 
aspect  ratio  and  Mach  number.  Using  the  airfoil  data,  we  estimate  the  section  lift-curve 
slope  as  Ci  =  6.0877/rad  for  a  NACA  653  -  018.  Mach  effects  are  introduced  based  on  the 
free-stream  Mach  number  for  the  calculation,  i.e.,  =  0.5.  The  maximum  thickness  of 

the  airfoil  occurs  at  x/c  =  40%,  which  gives  a  sweep  angle  of  Amax  t  =  20.7562°. 

The  MAC/C.G.  data  is  used  to  compute  the  moment  slope  via  Cm a  =  -CLahn-  Note 
that  the  center  of  gravity  is  estimated  using  relationships  between  the  center  of  gravity jmd 
the  location  of  aircraft  landing  gear.  Namely,  we  assume  the  landing  gear  is  located  15°  aft 
of  the  center  of  gravity.  The  finite-difference  values  in  the  tables  are  computed  using  the 
a  _  qo  ancj  a  —  g°  Euler  solutions.  An  argument  could  be  made  that  these  numbers  should 
be  computed  using  an  epsilon  change  in  the  angle  of  attack. 
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0.40 


Spanwise  Local  Moment  Coefficient 

YB-49  Flying  Wing,  M=0.5,  A=7.4,  1=0.25 


Figure  5.10:  Spanwise  moment  distribution  along  the  wing. 


CLa 

a 

Raymer 

Fin.  Diff.  SEM 

0° 

— 

—  4.3089 

4° 

4.2997 

5.0715  4.0818 

8° 

— 

—  4.0424 

Table  2:  Comparison  of  lift-curve  slope  between  semi-empirical,  finite-difference  and 
sensitivity-equation  method. 


5.1.9  Lateral/Directional  Stability 
Lateral/Directional  Stability 

Some  data  on  XB-35  lateral/directional  stability  is  presented  in  Northrop  in  his  Wright 
Memorial  Lecture  to  the  RAE.  For  this  pusher-propeller  configuration  a  considerable  amount 
of  weathercock  stability  is  provided  by  the  engine-drive  nacelles,  as  well  as  the  propellers 
themselves.  Since  our  model  does  not  include  these  features,  it  is  expected  that  the  bare- 
wing  estimate  for  Cn0  will  be  smaller  than  Northrop’s  data. 

The  sensitivity  of  the  wing  pressure  to  the  sideslip  angle  is  shown  in  Fig.  5.14  and 
Fig.  5.15.  The  pressure  sensitivity  on  the  upper  surface  of  the  right  wing  is  negative  leading 
to  an  addition  to  the  local  lift.  The  opposite  occurs  for  the  left  wing  where  pressure 
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Spanwise  Stability  Loading 
YB-49  Flying  Wing,  M=0.5,  A=7.4,  X=0.25 
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Figure  5.11:  Spanwise  stability  loading  distribution  along  the  wing. 


Cm* 

a 

Raymer  Fin.  Diff.  SEM 

0° 

_  _  -0.4613 

4° 

-0.4743  -0.6301  -0.4344 

8° 

_  _  -0.4450 

Table  3:  Comparison  of  pitching-moment  stability  between  semi-empirical,  finite-difference 
and  sensitivity-equation  method. 


sensitivities  are  positive  resulting  in  less  lift.  The  lift  differential  leads  to  rolling /yawing 
moment. 

A  comparison  between  Northrop  s  data,  finite  difference  and  the  sensitivity-equation 
method  is  given  in  Table  4  and  Table  5  for  the  rolling  and  yawing  rate  derivatives,  re¬ 
spectively.  These  are  sometimes  called  the  effective  dihedral  and  the  weathercock  stability 
derivatives.  The  trends  in  the  stability  derivatives  are  consistent  with  Northrop’s  data  as 
the  angle  of  attack  increases.  As  expected,  the  weathercock  stability  is  less  than  the  actual 
aircraft  which  uses  vertical  fins  for  increased  yaw  stability. 
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Figure  5.12:  Upper  surface 


Figure  5.13:  Lower  surface 

! 


5.1.10  Summary 

The  sensitivity  equation  method  (SEM),  applied  to  nonlinear  flow  analyses,  produces  a 
linear  boundary-value  problem.  The  solution  to  this  problem  is  a  flow  sensitivity  and  de¬ 
scribes,  in  linear  approximation,  the  dependence  of  the  flow-solution  on  a  scalar  parameter. 
The  derivative  of  any  non-linear  functional  of  the  flow  with  respect  to  this  parameter  can  be 
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Figure  5.14:  Upper  surface 


Figure  5.15:  Lower  surface 

1 


easily  computed  from  the  chain  rule.  This  latter  calculation  requires  evaluation  of  a  single 
inner  product  so  the  computational  effort  is  negligible. 

While  the  SEM  has  a  number  of  applications,  in  the  present  case  we  have  applied 
the  method  to  the  calculation  of  force  and  moment  stability  derivatives.  The  SEM  is 
computationally  efficient;  after  calculation  of  the  underlying  flow,  the  linear  boundary- 
value  problem  is  typically  solved  in  3%  of  the  time  required  for  a  non-linear  flow  solution. 
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~cT« 

a 

Northrop  Fin.  Diff.  SEM 

4° 

8° 

-0.029  -0.0347  -U.U318 
-0.057  -0.0696  -0.0667 

Table  4:  Computed  effective-dihedral  stability  derivatives. 


Cna 

a 

Northrop  Fin.  Diff.  SEM 

4° 

8° 

0.020  0.0159  0.0142 

0.026  0.0215  0.0203 

Table  5:  Computed  weathercock  stability  derivatives. 


A  single  sensitivity  can  be  used  to  calculate  the  derivative  of  a  number  of  force/moment 

^^'stabihty  derived ives^M  asymmetric  Bight  can  be  calculated  based  on  a  symmetric  flow 
solution'so^that1  configuration  symmetry  can  be  exploited.  The  method  naturaflyextends 
to  calculation  of  control  derivatives,  body-rate  derivatives  and  aeroelastic  derivatives. 

5.2  The  Impact  of  Finite  precision  Arithmetic  on  Sensitivity 

In  this  section  we  address  some  fundamental  issues  concerning  “time  marching”  numer- 

sidered  in  any  computational  study  of  non-unique  solutions  to  partial  diffiMOfcri 
that  govern  physical  systems  such  as  fluid  flows.  In  particular  numerical  calculations  have 
been  led  to  Suggest”  that  certain  Euler  equations  do  not  have  a  unique  solutiom  Fot 
Bmeers’  equation  on  a  finite  spatial  interval  with  Neumann  boundary  conditions  the  only 
steady  state  solutions  are  constant  (in  space)  functions.  Moreover 
ore.W  results,  for  any  initU d condition 

ut  numerical  steady  state  “solutions.”  These  erro¬ 
neous  solutiom  arise  out  of  the  necessary  flnite  floating  point  arithmetic  inherent  in i  every 
digital  computer  We  suggest  the  resulting  numerical  steady  state  solution  may  e 
^  a tlufcn  to  a  “nearby”  boundary  value  problem  with  high  sensitivity  to  chimges  m  the 
boundary  conditions.  Finally,  we  close  with  some  comments  on  the  relevance  of  this  ^tper 
to  some  recent  “numerical  based  proofs”  of  the  existence  of  non-unique  solutions 
equations  and  to  aerodynamic  design. 
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5.2.1  Introduction  and  Motivation 

During  the  past  twenty  years  there  have  been  a  significant  advances  in  computational  tools 
for  optimal  design  and  control  of  fluid  flows.  Many  of  these  tools  are  based  on  cascading 
simulation  software  into  optimization  and/or  control  schemes.  Although  this  method  has 
been  applied  to  some  complex  aerodynamic  control  problems,  very  little  has  been  done 
to  develop  a  rigorous  theory  of  convergence  for  the  resulting  optimal  design  and  control 
algorithms  Even  if  the  partial  differential  equation  that  defines  the  relationship  between 
state  and  control  variables  has  a  unique  solution  for  each  control,  numerical  approximations 
may  not  preserve  this  uniqueness.  In  particular,  the  discretized  state  equations  may  yield 
“false”  non-unique  solutions  which  may  drive  the  optimization  algorithm  to  an  incorrect 
design.  As  a  first  step  in  understanding  when  and  why  a  particular  algorithm  converges 
to  an  optimal  design  for  the  governing  partial  differential  equations,  it  is  essential  to  know 

when  and  why  these  extraneous  numerical  solutions  occur. 

Herewe  focus  on  just  such  an  issue.  In  particular,  we  use  a  simple  Neumann  boundary 
value  problem  for  the  one  dimensional  Burgers’  equation  to  illustrate  that,  because  of  finite 
precision  arithmetic,  a  convergent  numerical  algorithm  can  produce  false  (purely  numeri¬ 
cal)  solutions.  The  main  purpose  of  this  work  is  to  give  an  m-depth  examination  of  this 
model  problem  and  to  give  warning  in  the  use  of  numerical  based  proofs  of  uniqueness  for 

hydrodynamic  problems.  .  .  , . 

The  problem  discussed  in  this  paper  is  caused  by  computing  on  a  finite  precision  machine 

and  is  not  caused  by  roundoff  errors.  Also,  the  issue  considered  here  is  not  the  same  as 

supersensitivity  considered  by  other  authors.  _  , 

In  1993  Marrekchi,  while  working  on  a  Neumann  boundary  control  problem  for  Burgers 
equation,  observed  that  a  finite  element  scheme  used  to  design  feedback  controllers  produce 
non-constant  steady  state  solutions.  It  is  now  known  that  these  “solutions  are  purely 
numerical  and,  as  we  show  below,  most  numerical  methods  will  generate  such  false  solutions. 
If  a  boundary  value  problem  as  simple  as  Burgers’  equation  with  Neumann  conditions  can 
lead  to  such  complex  phenomena,  then  it  is  reasonable  to  expect  similar  difficulties  for 

potential  and  Euler  type  equations.  , 

Before  turning  to  Burgers’  equation,  we  present  a  simple  example  that  illustrates  the 

basic  difficulty. 

Example  1  Let  a  >  0  and  consider  the  initial  value  problem 

a  if  y(x)  <  1  +  10a 


^=s(y,a) 

ax  where  g{y,  a)  = 

1 2/(0)  =  1 

The  exact  solution  to  this  problem  is  given  by 

1  +  ax 

y{x)  = 


1  if  y(x)  >  1  +  10a 

,  0  <  x  <  10 

,  x  >  10 


1  4-  10a  +  (i  -  10) 

However  if  we  apply  a  standard  numerical  method,  such  as  Euler’s  method,  to  this  problem 
we  see  something  very  different  because  of  the  use  of  finite  precision  arithmetic.  Namely, 
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let  Ax  be  a  small  increment,  define  Xj  =  jAx  and  yj  ~  y(xj).  Then  the  Euler  iterates  are 


given  by 


yj+ 1  =  yj  +  g(yj  >  <*)  Ax,  yo  =  1 , 


and  for  j  =  0, 


yi  =  yo  +  y(yo,  <*)Ax  =  1  +  aAx. 


Now  for  aAx  sufficiently  small  and  with  finite  precision  arithmetic,  yi  =  yo  =  1,  be.,  if 
aAx  is  smaller  than  (l/2)2l~d  where  d  is  the  number  of  digits  (assuming  base  2  arithmetic), 
then  1 4-  aAx  =  1.  (Note:  the  smallest  positive  number  can  be  different  than  unit  rounding 
error  (machine  precision).)  Thus,  the  numerical  solution  to  this  problem  gives 


yj  —  1,  for  all  j  >  0. 

As  an  example,  on  a  particular  desktop  computer  using  MATLAB,  the  machine  pre¬ 
cision  is  given  in  a  variable  denoted  by  eps  whose  value,  for  this  machine,  is  eps  — 
2.220446049250313 10' 16  satisfying 


1  +  eps  =  1.0000000. 

As  a  numerical  experiment  set  Ax  =  .005  and  successively  solve  the  above  initial  value 
problem  using  Euler’s  method  as  described  above  with  a  sequence  of  decreasing  values 
a  _  j  =  1,  •  •  •  ,  14.  For  all  j  <  12  the  sequence  of  iterates  give  the  approximation  to 
the  correct  solution  but  when  j  =  14  the  iterates  give  yj  =  1.  We  note  that  in  this  case 

Ax  x  10~14  =  5  x  10-17  <  eps  <  Ax  x  10-12  =  5  x  10-15. 

The  solutions  for  these  cases  are  depicted  in  Figure  1.  The  lines  depict  numerical  solutions 
for  a  =  10-14  and  a  =  10-12. 


FIG.  1.  Machine  Floating  Point  Problem 


It  is  important  to  note  that  mesh  refinement  only  exacerbates  the  matter.  For  example, 
if  a  =  10-12  and  Ax  is  reduced  to  Ax  =  5  x  10-5,  then  Ax  x  10  12  =  5  x  10  <  eps. 

Therefore,  once  again  the  numerical  solution  yi  =  1  will  be  incorrect.  On  the  other  hand, 
for  a  =  0  the  solution  to 


! 


^=9(5.0) 

ax 


l  y( o)  =  i 


where  y(y,  0)  = 


if  y(x)  <  1 
if  y(x)  >  1 
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is  given  by  y(x)  =  1,  and  hence  the  numerical  solution  for  10-12  =  a  >  0  is  “correct"  for 
the  problem  defined  by  a  =  0.  We  suggest  that  some  of  the  numerical  problems  considered 
below  for  Burgers’  equation  follows  a  similar  pattern. 

Burgers’  equation  on  the  interval  (0, 1)  subject  to  Neumann  Boundary  Conditions  is 
given  by  the  dynamical  system 

wt  -  ewxx  +  wwx  =  0,  (5.19) 

x  G  (0, 1),  t  >  0,  e  >  0, 
wx(0,t)  —  wx(\,t)  —  0, 
w(x,0)  =  <j>(x). 

We  are  interested  in  the  corresponding  steady  state  problem 

-evxx(x)  +  v(x)vx(x)  =  0,  /5  2Qn 

vx(0)  =  u*(l)  =0. 

Solutions  of  (5.20)  are  called  stationary  (or  equilibrium)  solutions  of  the  unsteady  prob¬ 
lem  (5.19).  One  approach  to  the  development  of  numerical  methods  for  solving  (5.20)  is  to 
solve  the  time  dependent  problem  (5.19)  and  assume  that  *  v{-)  as  t  *  +oo  .  In 

order  to  construct  fast  and  accurate  “time  marching”  schemes  based  on  this  idea,  a  number 
of  points  must  be  considered.  In  particular,  one  should  address  the  following  issues. 

(a)  If  possible,  the  questions  of  existence  and  uniqueness  of  stationary  solutions  to  the 
boundary  value  problem  (5.20)  needs  to  be  answered.  These  are  still  open  questions 
for  many  fluid  and  gas  dynamic  problems. 

(b)  One  needs  to  know  that  for  reasonable  initial  data  </>(•),  the  time  varying  solution 
w(-,t)  exists  for  all  t  >  0,  lim  w{-,t)  =  v(-)  exists,  and  v(-)  is  a  stationary  solution. 

(c)  The  rate  at  which  w(-,t)  v(-)  is  important  because  it  can  influence  the  efficiency 

of  the  scheme. 

(d)  If  one  introduces  a  numerical  approximation  (with  spatial  mesh  size  Ax)  and  con¬ 
structs  the  numerical  solution  wAx(-,t)  with  the  property  that  as  Ax  — >  0  (i.e.  mesh 
refinement)  u/Ax(-,t)  — +  u>(-,£),  then  ^m^u;Ax(-,t)  =  uAx(-)  needs  to  exist. 

(e)  The  limit  uAx(-)  is  assumed  (or  proven)  to  be  a  good  approximation  to  u(-).  This 
issue  is  more  complicated  than  one  might  guess  and  it  can  fail  in  surprisingly  simple 
problems. 

Items  (a)  -  (e)  above  do  not  address  all  of  the  of  important  issues.  For  example,  as  we 
show  in  this  paper,  even  if  items  (a)-(d)  are  satisfied  and  one  can  prove  (this  means  using 
infinite  precision  arithmetic)  that  uAx(-)  -►  u(-),  then  problem  sensitivity  and  finite  precision 
arithmetic  can  produce  numerical  solutions  ^Ax(-)  that  do  not  approximate  any  stationary 
solution!  Thus,  it  is  possible  for  a  perfectly  sound  theoretical  algorithm  to  produce  “false” 
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numerical  solutions  to  the  steady  state  problem.  We  demonstrate  this  point  by  a  complete 
analysis  of  Burgers’  equation  with  Neumann  boundary  conditions. 

In  the  next  section  we  present  the  basic  analytical  results  for  systems  (5.19)  and  (5.20). 
We  then  present  some  numerical  results  to  illustrate  the  theoretical  results  and  to  demon¬ 
strate  a  numerical  anomaly  for  this  problem.  The  remainder  of  the  paper  is  devoted  to  the 
construction  and  analysis  of  numerical  schemes  so  that  items  (a)-(e)  are  established  and 
yet  numerical  solutions  generated  by  these  schemes  do  not  approximate  the  true  station¬ 
ary  solutions.  We  close  with  a  discussion  of  how  the  results  in  this  paper  may  relate  to 
“numerical  based  proofs”  of  nonunique  stationary  solutions  for  other  problems. 

5.2.2  Burgers’  Equation  with  Neumann  BCs 

As  pointed  out  in  our  earlier  work,  the  linearization  about  zero  of  (5.19)  in  L2( 0, 1)  is  the  one 
dimensional  heat  equation  with  Neumann  boundary  conditions.  A  well  known  consequence 
of  the  Fourier  representation  of  the  solution  for  this  problem  shows  that  the  unique  steady 
state  response,  for  any  initial  data  <f>  6  L?  (0, 1 ) ,  is  the  constant  function  with  value  equal 
to  the  average  of  the  initial  data,  i.e., 


For  Burgers’  equation  with  Neumann  boundary  conditions  it  is  easy  to  see  that,  w(x,  t )  = 
c  for  any  c  G  R  and  all  (x,t)  G  (0, 1)  x  [0,oo)  is  a  stationary  solution.  A  somewhat  deeper 
result  for  (5.19),  can  be  based  on  an  infinite  dimensional  version  of  the  Center  Manifold 
Theorem.  For  sufficiently  small  initial  data  in  Hx( 0,1 ),  the  solution  w(x,t)  of  (5.19)  tends 
to  a  constant  as  t  oo.  Since  the  Center  Manifold  Theorem  is  a  local  result  it  cannot  be 
used  to  make  any  general  statements  about  the  long  time  behavior  of  solutions  to  (5.19)  for 
larger  initial  data. 

Another  possible  approach  to  determining  the  long  time  behavior  of  solutions  to  (5.19) 
would  consist  in  determining  the  existence  and  properties  of  a  global  attractor.  Note  that 
since  such  an  attractor  must  contain  all  stationary  solutions,  it  must  contain  all  constants, 
so  it  must  be  unbounded. 

Actually,  a  more  relevant  first  question  to  answer  is  whether  solutions  even  exist  for  all 
time  for  alf  initial  data  <f>  G  L2(0, 1),  i.e.,  is  there  a  globally  defined  dynamical  system. 

In  the  1957  paper  by  Kiselev  and  Ladyzenskaya,  they  prove  the  global  existence  and 
regularity  for  multi-dimensional  Burgers’  equation  with  Dirichlet  boundary  conditions  using 
a  priori  estimates  and  the  maximum  principle.  While  this  paper  is  best  known  for  its 
contributions  to  Navier-Stokes  theory  it  also  has  a  section  devoted  to  Burgers’  equation. 

The  main  facts  needed  here  are  listed  in  the  following  jffieorem. 

Theorem  5.1  For  the  system  (5.19)  with  arbitrary  initial  data  <f>  G  L2(0,1)  and  0  <  T  < 
oo, 

a)  There  exists  a  unique  globally  defined  weak  solution  so  that  for  each  T  >  0 

w  G  C([0,T],  L2(0, 1))  n  L2([0,T],Hl(0, 1)), 
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b)  The  solution  is  instantly  infinitely  smooth  (and  therefore  classical)  for  t  >  0. 

c)  Therefore ,  there  is  a  globally  defined  dynamical  system  on  the  state  space  L2( 0, 1)  given 
in  terms  of  a  nonlinear  semigroup  {T*,£  >  0},  i.e.,  the  solution  is  given  by  w(x,t)  = 
Tt(<p)(x)  for  initial  data  <j).  This  semigroup  possesses  the  following  properties: 

i)  Tt  is  continuous  in  t  and  p  G  L2( 0, 1). 

ii)  Tt  is  compact  for  t  >  0. 

Hi)  There  exists  a  positive  continuous  monotone  increasing  function  a(£),  £  >  0 
such  that  a(0)  =  0  and 

|M<a(|M|),  <6(0,00),  V6L2(fi), 

which  means  that  the  system  is  globally  Lyapunov  stable, 
iv)  There  is  a  global ,  locally  compact  attractor  A. 

In  order  to  see  why  there  must  be  a  global,  locally  compact  attractor  we  note  that  by 
the  previous  theorem: 

1.  For  every  R  >  0,  the  ball  Ba(R)  is  an  absorbing  ball  for 

Br  =  G  L2{tt)  :  M  <  R} 

2.  This  implies  that  for  every  R  >  0,  the  dynamical  system  given  by  Tt,  restricted  to 
Ba(R)  has  a  nonempty,  compact,  connected  (local)  attractor,  which  we  denote  by  Ar- 

3.  It  is  clear  that  Ri  <  R2  implies  Ar,  C  Ar2,  and  hence  we  can  conclude  that  A  = 

U  Ar  is  the  global  attractor. 

R>  o 

4.  The  global  attractor  A  is  only  locally  compact  since  E  C  A. 

5.  Indeed,  for  R  >  0  sufficiently  small,  the  attractor  Ar  for  the  ball  Br( 0)  consists  only 
of  constants  (due  to  the  center  manifold  theorem),  i.e., 

AnBR(0)  =  AR  =  {c  :  |c|  <  a(R)}. 

We  should  comment  that  since  the  global  attractor  contains  all  stationary  solutions  and, 
as  we  have  already  mentioned  above,  every  scalar  is  a  stationary  solution,  the  attractor  is 
unbounded.  Due  to  Theorem  2.1  it  is  locally  compact.  The  exact  composition  of  the 
attractor  has  recently  been  settled  by  Edriss  Titi  and  Chdngsheng  Cao. 

Theorem  5.2  (Cao  and  Titi)  For  every  initial  cf>  G  L2(0, 1)  there  is  a  constant  c  so  that 

sup  \z(x,t)  —  c\~ — *  0. 

x6[0,l] 

The  dimension  of  the  global  attractor  is  one  and  consists  of  the  scalars. 
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For  fixed  e  and  for  small  initial  data,  numerical  approximation  of  the  solutions  to  (519) 
supports  the  conclusion  of  the  Center  Manifold  Theorem  and  Theorem  5.2.  In  particular, 
solutions  tend  to  a  constant  as  t  tends  to  infinity.  But  for  fixed  e  and  “certain”  initial 
data  (not  too  small),  some  numerical  solutions  converge  to  a  nonconstant  function.  These 
same  nonconstant  steady  state  limits  are  readily  obtained  using  many  different  numerical 
algorithms  and  on  various  computer  platforms.  We  are  led  to  conjecture  the  existence  of 
some  type  of  Numerical  Stationary  Solutions  for  the  problem  (5.19). 

One  class  of  initial  data  for  which  this  occurs  is  the  “antisymmetric”  functions,  that  is, 
functions  that  are  odd  about  x  =  1/2  in  the  interval  (0, 1), 

Las(°.  !)  =  {<?€  L2{ 0, 1)  :  <f>{x)  =  -<t>(  1  -  x)}  .  (5.21) 

For  initial  data  <j>  €  L2AS( 0,1),  a  straightforward  consequence  of  Theorem  2.1  is  that 
w(-,t)  G  L2as{0,  1)  for  all  t.  This  can  easily  be  seen  from  the  uniqueness  and  the  fact  that  if 
^  e  Lj lS(0, 1)  and  w(x,  t )  is  the  solution  of  (5.19),  then  the  function  z(x,t )  =  -w(l  -  x,t ) 
also  satisfies  (5.19)  and  hence 

w(x.t)  =  —w(  1  —  x,  t) 

i.e.,  w(',t)  €  L^s(0,  !)•  Note  that  a  continuous  function  (f)  in  L\ s(0, 1)  must  satisfy  </>(l/2)  = 
0  and  so,  for  t  >  0  a  solution  with  initial  data  <f>  €  L\ s(0,l)  will  satisfy  w{\/2 ,<)  =  0  for 
all  t>  0.  These  comments  establish  the  following  lemma. 

Lemma  5.3  The  Hilbert  space  L2AS{0, 1)  is  invariant  under  the  nonlinear  semigroup  Tt 
(defined  in  Theorem  5.1  part  d))  for  the  dynamical  system  (5.19).  Thus,  for  initial  data 
4>  ^  L2AS{Q,l), 

lim  w{x,  t )  =  0,  for  every  x  G  [0,  lj. 

£—•■00 

Note  that  a  weak  stationary  solution  must  satisfy  the  differential  equation 

(-«,,  +  y)  =0,  (5-22) 

where  all  the  derivatives  are  understood  as  weak  derivatives  and  the  equality  holds  for 
almost  every  x.  One  possibility  is  that  v  is  a  constant,  in  which  case  we  have, 


Clearly  a  constant  provides  a  stationary  solution  since,  in  addition,  it  satisfies  the  boundary 
conditions. 

The  only  distributional  solution  of  the  equation 

t//( x )  =  0 

is  a  constant,  so  any  other  stationary  solution  to  (5.19)  must  satisfy 

— evx  +  —  =  co,  coGM.  (5.24) 
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This  equation  can  be  solved  explicitly  with  solution  given  by 

v(x)  =  y/2co  tanh  (Cl  —  ’ 

where  c0  and  ci  are  arbitrary  constants.  A  straightforward  calculation  yields 

vx(x)  =  —~sech2  -  x))  > 


(5.25) 


(5.26) 


and  these  functions  cannot  vanish  at  x  =  0  or  x  =  1  (unless  Co  =  0).  Thus,  as  we 
already  know  from  Theorem  5.2  the  only  stationary  solutions  to  Burgers’  equation  satisfying 
homogeneous  Neumann  boundary  conditions  are  constants. 

The  nonconstant  functions  u(  )  given  by  (5.25)  form  a  two  parameter  family  depending 
on  the  parameters  co  and  cv  In  order  that  such  a  function  be  in  L\ s(0,l)  it  follows  that 
ci  =  1/2.  We  shall  focus  on  functions  h(-)  G  L\ s(0, 1)  defined  by 

h{x)  =  \/2co  tanh  (^(1/2  -  x))  .  (5-27) 


Although  h(-)  satisfies  equation  (5.22)  exactly,  it  only  approximately  satisfies  the 
ary  conditions  (to  within  exponentially  small  terms).  Namely,  the  functions  in  (5.27) 
(5.22)  and 


h'(x)  =  -^sech2 


(V2*> 
V  2e 


(1/2  -  xd  , 


bound- 

satisfy 

(5.28) 


which  for  small  e  and/or  large  co  gives 


ti{  0)  =  h'{  1)  = 


(5.29) 


where  o:  is  an  exponentially  small  positive  number. 

Thus,  if  a  is  close  to  0,  then  the  “nearby”  steady  state  problem 

equation 

— euXx(x)  +  v(x)vx(x)  =  0, 


defined  by  Burgers’ 
(5.30) 


with  (nonhomogeneous)  Neumann  boundary  conditions 

Ux(0)  =  Ux(l)  =  Cl) 


(5.31) 


will  have  non-constant  solutions  h{-)  given  by  (5.27).  As  we  see  below,  these  solutions  may 
appear  as  the  limit  (as  t  -*  +oc)  of  the  numerical  solutions  to  the  boundary  value  problem 

(5.19). 


5.2.3  Motivating  Numerical  Examples 

We  now  provide  several  examples  in  order  to  demonstrate  the  actual  behavior  of  numerical 
solutions  to  (5.19).  In  all  of  the  simulations  given  in  this  section  we  have  set  e  -  .1  and  have 
applied  the  finite  difference  method  presented  below  with  spatial  mesh  size  Ax  -  0.0125  - 
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1/80  and  temporal  mesh  size  At  =  0.0004.  We  conducted  these  numerical  experiments 
on  varying  time  intervals  and  for  a  variety  of  initial  data.  We  note  that  the  same  results 
occur  for  a  wide  variety  of  finite  element,  spectral  approximation,  and  other  finite  difference 
schemes. 

Observe  that  the  all  of  the  initial  data  </>(■)  belongs  to  L2AS{ 0,1)  so  that,  by  Lemma 
5.3,  the  solution  to  (5.19)  should  approach  zero  as  t  — ►  +oo.  Indeed,  this  is  exactly  what 
happens  when  the  initial  condition  </>(■)  is  “small”,  as  illustrated  in  Figure  1  and  Figure  2. 


FIG.  2.  Initial  data  <f>{x)  =  cos(nx)  and  final  time  T  =  10 


FIG.  3.  Initial  data  4>{x)  =  5(1/2  -  x)3  and  final  time  T  =  10 
However,  when  the  amplitude  of  the  initial  data  is  increased  we  observe  that  the  nu¬ 
merical  solutions  do  not  converge  to  zero.  In  fact,  they  appear  to  converge  to  a  solution  of 
the  nearby  problem  defined  by  (5.30)-(5.31).  More  will  be  said  about  this  phenomenon  in 
the  conclusions.  In  Figure  4  we  have  increased  the  amplitude  of  the  cosine  function  from  1 
to  5  and  in  Figure  5  we  have  increased  the  amplitude  of  the  cubic  initial  function  from  5  to 

20 
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FIG.  5.  Initial  data  <f>(x)  =  50(1/2  -  x)3  and  final  time  T  -  0.25 


Next  we  consider  an  initial  condition  to  show  that  the  same  results  hold  even  for  initial 
data  that  are  not  strictly  decreasing.  In  particular,  we  consider  the  initial  function  4>{x )  - 

v  /  ^  V 


A  (  -  -  x 
K  4 


- x 


-x]  for  various  values  of  A. 


FIG.  6.  Initial  data  <j>(x)  =  5 
final  time  T  =  10 


(M(M0Hand 


FIG.  7.  Initial  data  <f>(x)  =  50  -  x^  x 


and 

final  time  T  =  0.25 

Once  again  we  see  that  for  A  small,  as  t  -»  oo  the  numerical  solution  converges  to  zero, 
as  it  should.  But  for  a  larger  amplitude  A  the  numerical  solution  converges  to  a  nonconstant 

steady  state. 


34 


Finally  we  take  the  negative  of  the  initial  data  that  gave  rise  to  numerical  solutions 
which  converged  to  nonconstant  steady  states  and  we  note  that  these  solutions  ten  o 
zero  very  quickly.  Consequently,  nonconstant  numerical  steady  state  solutions  do  not  occur 
simply  because  of  larger  magnitude  initial  data. 


X 


FIG.  8.  Initial  data  </>(x)  =  -5cos(7rx)  and  final  time  T  -  2.5 


X 

FIG.  9.  Initial  data  <j)(x)  =  -50(1/2  -  x)3  and  final  time  T  =  2.5 


1  t 


FIG.  10.  Initial  data  <£(x)  =  — 
and 

final  time  T  =  2.5 
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initial  data 
cos(7r(l/2  —  x)) 
5cos(7r(l/2  —  x)) 
— 5cos(7r(l/2  —  x)) 
5(1/2 -x)3 
50(1/2  -  x)3 
-50(1/2 -x)3 


T 

co 

10.00 

5.85e-06 

0.50 

1.24e+01 

2.50 

8.71e-04 

10.00 

1.74e-07 

0.50 

1.56e+01 

2.50 

5.44e-04 

10.00 

7.54e-08 

0.50 

7.31e+00 

2.50 

2.84e-04 

||/l-£||oo 

h'(  0) 

8.57e-04 

-5.85e-05 

8.08e-02 

-7.65e-09 

5.32e-02 

-8.62e-03 

2.74e-05 

-1.74e-06 

1.15e-01 

-4.64e-10 

4.29e-02 

-5.41e-03 

1.18e-05 

-7.54e-07 

3.56e-02 

.  -1.46e-06 

3.26e-02 

-2.83e-03 

TABLE  I:  Initial  conditions,  final  time/T,  numerically  computed  co,  maximum 
difference  between  h(-)  and  hj  =  h(jAx)  for  j  =  0, 1,  •  •  •  ,  (N  4-  !)• 


In  order  to  draw  attention  to  the  actual  functional  values  of  the  apparent  nonconstant 
steady  state,  in  the  above  examples  we  denote  by  h(x)  the  numerically  computed  stationary 
solution.  To  be  precise, 

h(x)  =  lim  wAx(x,t ), 

t— KX> 


where  wAx{x)t) 
We  now  let 


is  computed  by  the  Crank-Nicolson  scheme  described  in  the  next  section. 


co  = 


-ehx{0)  + 


#(0) 

2 


In  every  Figure  2-10  the  curves  (plotted  on  the  right  side)  depicting  the  numerical  solu¬ 
tion  at  the  final  time  value  T  actually  contain  plots  of  both  the  numerical  solution  />(•)  and 
the  explicit  function  defined  by  equation  (5.27)  which  only  “nearly”  satisfies  the  boundary 
conditions.i.e.,  ti(0)  =  h'{  1)  =  -a.  Observe  that  h(-)  and  h(-)  are  not  distinguishable  on 
this  scale.  In  Table  I  we  have  compiled  a  list  for  comparison  of  various  parameters  for  our 
numerical  examples.  For  all  examples  we  have  set  e  =  .1  and  N  =  80  (the  number  of  nodes 
for  the  Crank-Nicolson  scheme). 


5.2.4  A  Symmetrized  Crank-Nicolson  Method 

It  has  been  observed  that  applying  standard  numerical  methods  to  the  system  (5  19)  for 
certain  initial  data  in  L^s(0, 1)  produces  numerical  nonconstant  steady  state  solutions. 
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However,  Lemma  5.3  implies  that  not  only  must  these  solutions  be  constant,  the  constant 
must  be  zero.  The  objective  of  this  section  is  to  prove  that  the  observed  difficulty  is  due 
to  the  use  of  finite  precision  arithmetic.  To  do  this  we  need  to  first  obtain  a  very  stable 
numerical  scheme  on  which  to  base  our  proofs.  The  reason  for  this  concern  is  that  if,  for 
example  we  employ  a  straightforward  finite  difference  scheme  to  approximate  (5.19)  on  the 
whole  interval  (0,1),  then  no  matter  how  small  the  spatial  discretization  or  the  time  step 
size,  eventually  round-off  error  will  corrupt  the  data  at  x  =  1/2  (where  the  solution  is 
known  to  be  zero  for  all  time).  Because  of  this,  there  will  always  be  a  time  at  which  the 
numerical  method  will  undergo  a  rapid  change  and  then  generally  converge  to  a  nonzero 
constant.  The  sign  of  this  constant  depends  on  whether  the  value  of  the  numerical  solution 
first  begins  to  drift  positive  or  negative  at  x  =  1/2.  So  our  first  step  in  obtaining  a  more 
stable  numerical  scheme  consists  of  the  reduction  to  a  problem  with  a  Dirichlet  boundary 
condition  at  x  =  1/2. 

Prom  Lemma  5.3  we  see  that  for  initial  data  <j>  €  L2AS(0, 1)  we  can  replace  problem  (5.19) 
by  the  system 

wt  —  ewxx  +  wwx  =  0,  (5.32) 

*€(0,1/2),  t>  0, 

Wi(0,t)=0,  w(l/2,t)  =  0, 
w(x,0 )  —  <p(x). 


Thus,  we  need  only  solve  the  Burgers’  equation  on  an  interval  of  half  the  length  and,  more 
importantly  for  numerical  calculations,  we  have  replaced  a  Neumann  boundary  condition 
at  x  =  1  by  a  Dirichlet  Boundary  condition  at  x  =  1/2. 

Consider  a  standard  implicit  finite  difference  scheme  for  the  system  (5.32)  on  the  interval 
(0,1/2)  with 

X{  —  zAx,  i  —  0,  1,  2,  ...  ,  N,  Ax  —  -rjr , 


A  t 

7  = 


2Ax' 


K 


At 

(Ax)2 


,  tj  =  jAt,  j  =  0,1,2, 


and,  in  a  standard  notation,  define  the  approximation 


w 


{xi,tj),  i  =  0, 1,2, . . .  ,N,  j  =  0,1,2, - 


Then,  we  have  for  i  =  1, 2,  •  •  •  ,  (N  —  2) 

wi,j+l  —wij  +  2  K+l-i  —  ^wiJ  T 

fx  K 

+  ^  [wi+ij+i  -  2wi,j+i  +  ^i-ij+i]  (5.33) 

+  -^WiJ  \wi+lj  -  Wt-lj]  +  2Wij+l  lu,i+1d  _  > 
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(5.34) 


and  for  i  =  0  and  i  =  N  —  1  we  have 

woj+i  =woj  +  —  [2wij  —  2^oj]  +  2  [2^1  ,j+i  —  2u,oj+i]  i 
lOjV-lj+1  =WA-l,j  +  -  [-2WN-1J  +  WN-2,j] 

+  ^  [— 2u>jv-ij+i  +  WN-2,j+i]  (5-35) 

7  7 

+  -WN-ijWN-2,j  +  — 2,j+l* 

Note  that  to  obtain  (5.34),  we  let 

(Wl  j-w_hj)/( 2Ax)  =  0  and  uw,j  =  0  for  all  j  >  0. 

For  this  difference  scheme  let 

W-N+l,j 
w—N+2,j 
Wj  =  ; 

™N-2,j 
.  WN-lj  . 

denote  the  solution  at  time  step  tj.  Then,  we  shall  prove  the  following  result. 

Theorem  5.4  For  any  piecewise  continuous  initial  data  w(x,0)  €  L2AS(0,l)  and  N  suffi¬ 
ciently  large 

IKIh^o. 

That  is,  for  every  initial  condition  in  L2AS( 0, 1),  the  finite  difference  solution  converges,  for 
every  N  sufficiently  large,  and  the  resulting  limit  is  zero. 

For  analysis  purposes,  it  is  useful  to  extend  this  difference  scheme  to  the  symmetric 
interval  (-1/2, 1/2)  with  Dirichlet  conditions  at  each  end  as  long  as  we  note  that  in  the 
approximations  of  the  first  derivative  terms  the  signs  must  be  reversed.  Thus,  for  i  = 
-N  +  2,  -N  +  3, . . . ,  -1  we  have 

Wij+i  —Wij  +  —  [wi+i,j  —  2  Wij  +  Wi-ij] 

+  ^  ,j+i  —  +  W{— lj+i]  (5.36) 

-  ^Wij  [wi+lj  -  Wi-lj]  -  7£wi,j+ 1  [Wi+l,j  -  > 

4 

with 

w-ZV+ij-f i  —w-N+ij  +  2  [”2w_tv+ij  +  iv ~N +2,j] 

+  £  [-2W-N+IJ+1  +  U>-JV+2j+l]  (5-37) 

7  7 

-  ±W-N+l,jW-N+2,j  -  - W-N+hj+lw-N+2,j+l- 


38 


The  resulting  difference  scheme  can  be  written  in  vector  form  as 

Awj+i  =  Bwj  +  |w j  +  ^Wi+ 1, 


where 


WJ  = 


Wj+1  = 

After  simplification  this  system  can  be  written  as 


W-M+l,jw-N+2J 
W-N+2J  (W-N+l ,j  -  W-N+3,j) 

WN-2,i  {wN-3,j  -  wN-l,j) 
U)N-l,jWN-2j 

W-M+l,j+lw-N+2J+l 
W-N+2J+1  {w-N+l,j+l  ~  W-N+3j+l) 


WN-2J+1  (WN-3J+1  ~  WjV-lj'+l) 
WN-l,j+lwN-2,j+l 


Awj+1  =  -Bwj  -  lWj  -  JWj+1, 


7, 


with 


A  = 


-2-  - 
K 


— 2 - 1 

K 


B  = 


2-2 

K 

1 


0 


(5.38) 


1  "t"  K 

K 

_  2 

K 

~2 

1  4*  ft 

K 

~2 

0 

,  B  = 

1  -  K 

K 

2 

K 

2 

1  -  K 

0 

K 

2 

0 

K 

~2 

1  +  K 

0 

K  1 

2  '“"J 

(5.39) 


- -2  1 
K 


1  --2 
K 


(5.40) 


0  1  -2  - 

It  is  well  known  that  the  matrix  A  is  invertible  and  so  once  again  we  can  rewrite  the 
system  (5.40)  as  _ 

w,+1  =  -A~lBwj  —  — A_1Wj  —  — A-1Wj+1)  (5.41) 

and  thus  we  have 

||wJ+i||2  <  || -A- 1 II2 !l wy II2  +  -||^_1||2||Wj||2  +  ^-\\A  1  IbllWj+xlU-  (5-42) 

We  need  the  following  result  in  order  to  prove  Theorem  5.4. 
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Lemma  5.5  The  eigenvalues  of  the  (m  -  1)  x  (m  -  1)  tri-diagonal  matrix 

'a  1  O' 

1  O'  1 


0  1  a] 


are 

(2  +  a)  -  4  sin2 

Thus,  the  eigenvalues  of  A  are 


i  =  1,2, ...  ,2m  -  1. 


i  =  1,2, . . .  ,2N  —  1, 


and  the  eigenvalues  of  A  1 J5  are 


2  A  .  2 
- 4  sin 

K 


Z7T 

4N 


2  ,  .  2 
- 4  sm 

K 


m 

AN 


i  =  1,2, . . .  ,2iV  —  1. 


Proof  (of  Theorem  54)  Let 

IlyT1  Silo  =  8  =  max  fa,  and  notice  that  /?  <  1. 
11  ^  1<K2JV-1 


Since 


we  have 


IIA-'lb  < 

II  Wj+i  II2  =  /3||wj||2  +  |l|Wj||2  +  |||Wj+i||2, 


and  the  inequalities 

l|Wj||2  <  2 llwjlil,  ||Wj+1||2  <  2||wj+i||l, 


(5.43) 


(5.44) 


imply  9 

IK+i||2  <  /?||wj||a  +  7IKII1  +  7l|wi+i||5. 

If  C  =  l|wJ+1||2,  s  =  ||wj ||2 >  then  (5.45)  can  be  written  as 

C  <  /?s  +  7s2  +7C2. 

Define  r  =  /3s  +  7s2  and  so  7C2  —  C  +  r  >  0  so  that  C  must  satisfy 


(5.45) 


C< 


1  ~  \/l  ~  47^  <  (47r/2)  +  (2/9) (47r)2  for  4~r  <  3^ 


27 


27 


(5.46) 
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In  order  to  justify  the  last  inequality,  we  show  that  for  0  <  |  <  3/4 

t  ,  2£2 


This  inequality  follows  from  the  chain  of  implications 

3 


0<£< 


1  +  M 

2  9 


*G)' 


Rearranging  and  multiplying  by  £2,  gives 

-^2  +  (f  +  ^2)  so' 

Adding  (1  -  £)  to  both  sides  and  rearranging  terms  gives 

The  left  side  of  this  expression  is  a  perfect  square  so  we  have 


and  finally  we  have 


4  ,  2£2 


1-VT^<!+  9 


Setting  C  =  4 qr  yields  inequality  (5.46). 

Returning  to  the  estimate  (5.46)  we  have  that  if  4qr  <  3/4,  then 

167  2 

C  <  r  +  -gV. 

Now  choose  d  so  that  (0  +  d)  <  1.  Suppose  that,  as  a  very  conservative  estimate, 
s  =  || Wj H2  <  5/(47).  Then 


r  =  0s  +  7  s' 


*H) 


and 


C<  £+7  *  + 


I67 


=  {i3+a)s  + 

<  ( 0  +  d )  s. 


a  d\  a 
p  +  l) s ^ 

S 


s8H) 
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||w0||2  < 


Hence,  if  ||wj||2  <  5/(47)  then 

||wj+i||2  <  (^  +  5)||wj||2,  where  /3  +  5<l. 

This  implies  that  ||Wj||2  — ►  0  as  j  — ►  oo,  provided  that  ||w0||2  <  5/(47).  Note  also  that 
when  || Wj ||2  <  5/(47),  then  47r  <  d(0  +  5/4)  <  5.  Thus,  we  also  need  5  <  3/4. 

Recalling  the  definition  for  7  =  At/ (2 Ax)  we  have 

„  „  ^  5  ^  (l-/3)Ax 

S  4?  5  2Af~ 

Generally  we  would  require  that  At  is  somewhat  smaller  than  Ax.  Indeed,  it  is  possible 

‘°take  ..  e(Ax)2  ^  e(Ax)2 

^  4  <  2 

so  that  the  initial  requirement  would  be 

2(1  -£)=  4(1  -  0)N_ 
eAx  e 

where  in  the  last  equality  we  have  used  Ax  =  ^.  Therefore,  we  see  that  this  condition 
is  not  really  a  restriction,  except  for  numerical  difficulties,  since  we  can  take  N  sufficiently 
large  to  include  any  given  piecewise  continuous  initial  data  and  this  completes  the  proof.  0 

5.2.5  Stationary  Solutions  of  Explicit  Finite  Difference  Schemes 

In  this  section  we  consider  several  explicit  finite  difference  schemes  for  the  problem  (5.19) 
and  make  some  observations  concerning  the  existence  of  nonconstant  stationary  solutions 
for  these  discrete  schemes.  In  this  section  we  assume  that  N ,  the  number  of  spatial  nodes, 
is  even  and  we  set  Ax  =  1  /N,  for  i  =  0,1,...,  Xj  =  iAx  and  tj  =  j  At.  We  consider 
three  different  schemes  based  on  forward  difference  in  time,  central  difference  for  the  second 
order  term,  and  three  different  approximations  for  the  convective  term.  (Recall  the  notation 
Wij  «  w(Xi,tj).) 

Example  2  Centered  Difference  wwx :  For  i  =  0, 1,2,  •  •  •  ,N 

(5-47) 


{Wij+l  -  Wij ) 

i,j  ij  Wi— ij 

-  WU 

Wi+lj  -  Wi-lj 

- - — -  =  € 

At 

(Ax)2 

(2Ax) 

IW-IJ  =  Wij  WM+Uj  =  WN-l  ,j- 

The  steady  form  of  (5.47)  is  given  in  terms  of  a  function  v(x)  with  V{  =  vix/)  by 

vi+i  -  2 Vi  +  Vi- 1  -  riVi  (vj+i  -  Vi-i)  =0,  i  =  0,...,N 

Ax 

v-i  =  vi,  Vjv+i  =  vn-1  and  rq  =  — . 


(5.48) 

(5.49) 


For  this  case  (5.48)  has  one  stationary  solution,  the  constant  solution,  i.e.,  for  any  c  £  M, 
Vi  =  c  for  i  =  0, 1, . . . ,  N  is  a  solution.  But,  there  is  another  solution  given  by 

'0,  i  =  N/2 

Vi  =  <l/ri,  0  <  i  <  N/2  —  1  (5-50) 

-1/n,  N/2  +  1  <  i  <  N. 
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Example  3  Centered  Difference  (w2)x/2:  For  i  =  0,1,2,--  ,  N 


(Wij+l-Wij)  " 

2  “f“  Wi—lJ 

l 

wi+ 1  j  ”  wi-lj 

At 

1 - 

> 

to 

K _ 

~  2 

(2Ax) 

W-l,j  =  WXj  WN+lj  =  WN-1J- 


(5.51) 


The  steady  form  of  (5.51)  is 


vi+ 1  -  2vi  +  Vi-i  -  r2  (vt2+1  -  =0,  i  =  0,...,N 

Ax 

V_1  =  Vi,  VN+ 1  =  VN- 1  and  r2  =  — . 


(5.52) 

(5.53) 


For  this  case  (5.52)  has  a  solution  Vi  =  c  for  i  =  0, 1, . . . ,  N  for  any  c  €  R.  Again  there 
is  a  second  solution  given  by 


0, 

1/^2, 

-l/ra, 


i  =  JV/2 

0  <  i  <  N/2  -  1 
JV/2  +  1  <  *  <  AT. 


(5.54) 


Example  4  Forward  Difference  wwx:  For  i  =  1, 2, ,  N  —  1 


(Wij+l  -  Wij)  _  " 

2  Wij  +  Wi—ij 

~  Wij 

Wij  -  Wi-lj' 

At 

(Ax)2 

(&x) 

looj  =  wij  -  wx-ij- 


(5.55) 


The  steady  form  of  (5.55)  is 

Vj+l  -  2vi  +  Ui_i  -  r3Vi  (vi  -Ui_i)  =  0,  *  =  (5.56) 

v0  -  vi,  vx  =  vx-i  and  r%  =  — .  (5.57) 

However,  in  this  case,  one  can  show  that  the  only  stationary  solution  is  the  constant  solution 
Vi  =  c  for  i  =  1 , . . . ,  N  —  1 . 


The  above  examples  illustrate  that  discrete  versions  of  the  steady  state  problems  can 
have  non-constant  discrete  solutions.  Thus,  if  one  uses  one  of  the  methods  (5.47)  or  (5.51), 
then  any  numerical  algorithm  (time  marching,  direct,  etc.)  could  produce  a  discrete  non¬ 
constant  stationary  solution.  This  can  happen  even  when  the  original  partial  differential 
equation  does  not  have  such  nonconstant  stationary  solutions.  Thus,  the  numerical  station¬ 
ary  solutions  are  not  approximate  solutions  to  the  steady  state  partial  differential  equation 
(plus  boundary  conditions). 

On  the  other  hand,  the  scheme  (5.56)  only  has  constant  discrete  stationary  solutions. 
Therefore,  one  might  expect  that  if  the  discrete  equations  (5.56)  are  used  then  numerical 
solutions  based  on  this  type  of  algorithm  will  not  produce  false  solutions.  However,  because 
of  finite  precision  arithmetic,  this  assumption  is  not  valid. 


43 


5.2.6  The  Effect  of  Finite  Precision  Arithmetic 

In  reality,  all  calculations  are  done  in  finite  precision  arithmetic.  Even  calculations  using 
computer  algebra  systems  that  purport  to  be  capable  of  infinite  precision  are  actually  limited 
by  memory  and  storage  limitations.  More  realistically,  floating  point  arithmetic  is  commonly 
used  in  computational  work  and  floating  point  arithmetic  is  based  on  a  finite  set  of  numbers 
and  a  finite  precision  arithmetic.  Furthermore  different  machines  have  a  different  set  of 
numbers  and  precision.  We  plan  to  show,  by  way  of  examples,  that  the  reason  for  the 
anomaly  observed  in  this  work  is  due  to  the  necessary  use  of  finite  precision  arithmetic.  We 
have  already  seen  in  the  last  section  that  using  exact  arithmetic  the  symmetrized  Crank- 
Nicolson  numerical  scheme  must  converge  to  zero  for  initial  data  in  L\ s(0, 1).  In  this  section 
we  show  that  by  altering  only  the  magnitude  of  the  initial  data,  the  value  of  the  viscosity 
e  and  the  precision,  we  can  generate  solutions  that  converge  either  to  zero  or  to  one  of  the 
analytic  solutions  that  only  approximately  satisfy  the  boundary  conditions.  We  note  that 
for  e  fairly  large  (for  example  e  =  1/2)  convergence  to  a  nonconstant  numerical  stationary 
solution  requires  larger  magnitude  initial  data  and  for  smaller  e  we  can  take  the  magnitude 
to  be  much  smaller.  (Note  that  the  Crank-Nicolson  scheme  is  used  for  all  the  calculations 
in  this  section.) 

In  our  first  example  we  consider  an  initial  condition  which  happens  also  to  correspond 
to  a  possible  fixed  point  of  the  finite  centered  difference  scheme  used  in  Example  2.  Fix  N, 
e,  and  U  >  0,  then  define  the  initial  data 

fu ;  o  <  x  <  1/2 

cj)(x)  =  <  0,  x  =  1/2  .  (5.58) 

1/2  <  x  <1 

Although  we  do  not  include  the  proof  here,  we  have  shown  that,  with  full  (infinite)  preci¬ 
sion  arithmetic,  the  spatially  centered  difference  spatial  discretization  and  Euler  marching 
scheme  for  U  <  l/r\  =  2iVe  converges  to  zero  with  increasing  time.  Thus,  this  example 
provides  a  good  test  for  our  hypotheses  that  numerical  stationary  solutions  arise  from  finite 
precision  arithmetic. 

In  this  numerical  example,  we  take  N  =  40  and  e  =  1/5  so  that  in  (5.50)  r\  =  1/16. 
Thus,  a  stationary  solution  for  the  centered  difference  method  is  given  by 

'l/n,  0  <  x  <  1/2 

v(x)  =  <  0,  x  =  1/2  .  (5.59) 

^ — 1/ri,  1/2  <  x  <  1 

For  the  first  numerical  run  we  take  initial  data  with  U'  =  8  and  we  successively  choose 
the  number  of  significant  digits  to  be  d  =  4, 8, 16  for  numerical  simulation. 
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X 


d=16 

FIG.  11.  Piecewise  Constant  Initial  Condition  U  =  8 


Notice,  as  shown  in  the  following  table,  that  the  values  continue  to  drop  for  d  —  16 
while  a  steady  state  is  reached  for  d  —  4  and  d  =  8. 


t  d  =4 

0  8.00000000000000 
10  7.99900000000000 

20  7.99900000000000 

30  7.99900000000000 

40  7.99900000000000 

50  7.99900000000000 

60  7.99900000000000 

70  7.99900000000000 

80  7.99900000000000 

90  7.99900000000000 

100  7.99900000000000 


d  =8 

8.00000000000000 

8.00000000000000 

8.00000000000000 

8.00000000000000 

8.00000000000000 

8.00000000000000 

8.00000000000000 

8.00000000000000 

8.00000000000000 

8.00000000000000 

8.00000000000000 


d  =16 

8.00000000000000 

7.99985129822842 

7.99970174625331 

7.99955214983113 

7.99940250894646 

7.99925282355433 

7.99910309364852 

7.99895331919820 

7.99880350016796 

7.99865363654645 

7.99850372828631 


TABLE  II:  Values  at  x  =  0 


Now  consider  a  smaller  value  [7  =  4  and  repeat  the  same  calculations. 
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0  0.2  0.4  0.6  0.8  1 


d  =  16 

FIG.  12.  Piecewise  Constant  Initial  Condition  U  =  4 


to 


Notice  as  shown  in  the  Figure  11  and  the  following  table,  the  computed  values  continue 
decrease  to  zero  for  d  =  8  and  d  =  16,  but  a  non-constant  solution  is  reached  for  d  =  4. 


t  d  =4 

0  4.00000000000000 

10  4.00000000000000 
20  4.00000000000000 
30  4.00000000000000 

40  4.00000000000000 

50  4.00000000000000 

60  4.00000000000000 

70  4.00000000000000 

80  4.00000000000000 

90  4.00000000000000 

100  4.00000000000000 


d  =8 

4.00000000000000 

3.58872030000000 

2.13481320000000 

0.00000025606987 

0.00000000000000 

0.00000000000000 

0.00000000000000 

0.00000000000000 

0.00000000000000 

0.00000000000000 

0.00000000000000 


d  =16 

4.00000000000000 

3.58892931217814 

2.13685522204510 

0.00000025848659 

0.00000000000000 

0.00000000000000 

0.00000000000000 

0.00000000000000 

0.00000000000000 

0.00000000000000 

0.00000000000000 


TABLE  III:  Values  at  x  =  0 

If  we  now  take  U  =  1/r  =  16  then  the  initial  condition  is  a  stationary  solution  for  the 
particular  values  N  =  40  and  e  =  1/5.  Indeed,  when  we  run  the  same  calculation  as  above 
we  obtain  the  plot  given  in  Figure  13. 
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FIG.  13.  Relation  to  Stationary  Solution  of  Discrete  Problem 

In  each  case  above,  at  each  time  step,  we  have  also  computed  the  numerical  value  of  the 
derivative  of  the  solution  w  at  x  =  0  and  used  the  fact  that 

w(0,t)2  wx(0  ,t) 

co  -  g  C  2 


to  obtain 


,m2 


Co 


w. 


oj 


(  W±j_-_WQj\ 

{  Ax  )■ 


At  each  time  step  the  above  graphs  also  include  the  graph  of 


whose  derivative  at  zero  is 

p 

wx(0.  tj)  «  t^jx(O)  =  — ^-sech 

In  the  last  case  we  have  cj  =  128.0  for  all  j  and  the  approximation  to  the  derivative  is 

WjX(0)  =  — 1.09 10-14 


5.2.7  Finite  Arithmetic  and  Convergence  to  Steady  State 

Let  0  denote  the  base  for  a  computer  system  and  t  the  number  of  digits.  On  the  interval 
[0m~l,0m],  the  floating  point  numbers  are  evenly  spaced  with  separation 

pm-1  a,  «2  er 

I - 1 — 1 - 1 

FIG.  14.  Floating  Point  Numbers  =  oli  +  0m~t 
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Let  x\,X2  G  [/3m  1 ,0m]  be  two  floating  point  numbers.  If  |xi  —  I2I  <  7,^™  *>  t^ien 
xi  =  X2-  Thus,  if 

or,  if 

\xi  ~xl\  ^  1  fli-t 
1*1 1  2**  ’ 

then  xi  —  X2- 

Consider  the  non-homogeneous  problem  for  Burgers’  equation  on  the  interval  [0, 1/2] 
given  by 

wt  —  ewxx  4-  wwx  =  0,  (5.60) 

wx(0,t)  =  -a,  a  >  0, 
tw(l/2,t)  =0 
w(x,0)  =  4>(x), 


If  a  =  0  then  we  know  from  Theorem  5.2  that  w(x,t)  — ♦  0  as  t  — ►  00.  If  a  7^  0  we  expect 
from  our  numerical  evidence  (and  the  discussion  in  the  next  section)  that 

w(x,  t )  h(x)  =  \/2co  tanh  ^^^(1/2  —  x)^  , 

h'(0)  -  jsech2  =  -«• 


where 


Consider  the  Crank-Nicolson  scheme  for  (5.60)  and  let 

(WU  -  w-i j)  _ 

2Ax 


so  that 

It  follows  that 


wij  =  — 2a(Ax) 


2aAx 

2aAx' 

w-lj  =  WlJ 

1  + 

~  w\  J 

1  + - 

WU  J 

Wqj  _ 

and  hence,  if  2a(Ax)/w0J  <  (l/2)/?1_t,  then  w-ij  =  wij.  This  implies  that  the  condition 

Kj  ~  w-i j)  _  n  . 

2Ax 


is  equivalent  to 


if 


(u-i,j  -  w-u) 

2Ax 


2a  (Ax) 
wo  j 
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Thus,  even  when  a  -  0  the  solution  on  any  computer  may  converge  to  (approximately) 
h(x). 

In  order  to  examine  this  more  closely,  consider  the  inequality 

2<*(Ax)  l-i  -t 

m,j  2 

It  is  reasonable  to  expect  that  v.'o j  ~  0(0)  if  a  nonzero  steady  state  is  quickly  reached. 
Consequently,  if 

^>>1, 
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then 


(f)(0)  «  h( 0)  «  \/2 cq. 


fey  — 

Assuming  that  ^  0  1,  we  have 
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a  «  4 —  exp 
6 


26 


and 


2a(Ax)  <  l^i_t 
wo  j  2 


can  be  approximated  by 


8(Ax)tf>(0) 


c 


,  *0)\ 

2e  eXPl^J 


<  2^"- 


(5.61) 


Consider  once  again  the  numerical  example  above  with  initial  data  given  by  (5.58) 
restricted  to  the  interval  [0, 1/2].  In  particular,  we  take 


4>(x)  = 


'U,  0  <  x  <  1/2 
10,  x  =  1/2 


(5.62) 


We  take  d>(0)  =  U,  the  number  of  digits  d  =  t,  f3  =  10,  5x  =  1/80,  e  =  1/5  and  for  the 
numerical  scheme  we  have  taken  a  time  step  size  of  6t  —  1/4000.  If 


U  f-5U\  l1ftl_d 

4eXPl'2“j<21°  ' 


we  expect  the  computed  solution  may  approach  a  nonzero  steady  state. 

For  U  =  4  and  d  =  5  we  have 

^exp  (-rr)  =  0.4539992976248510-4  <  0.5  x  10“4  =  ^lO1-'*, 
so  we  expect  that  the  solution  will  converge  to  a  nonzero  steady  state. 
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FIG.  15.  Trajectories  at  tj  =  10 j,  j  =  0, 1,  ■  •  •  ,10  with  d  =  5 
For  (7  =  4  and  d  —  6  we  have 

^  =  0.4539992976248510-4  >  0.05  x  10-4  =  ^101_d 

so  we  expect  that  the  solution  will  converge  to  a  zero.  This  behavior  is  demonstrated  in 
Figures  15,  16  and  Table  IV. 


FIG.  16.  Trajectories  at  tj  =  10 j,  j  =  0, 1,  •  •  •  ,  10  with  d  =  6 


t 

d  =  5 

10.00 

4.000000 

20.00 

4.00000 

30.00 

4.00000 

40.00 

4.00000 

50.00 

4.00000 

60.00 

4.00000 

70.00 

4.00000 

80.00 

4.00000 

90.00 

4.00000 

100.00 

4.0000 

d  =6 

3.60426000000 

2.32885000000 

0.000000737544 

0.000000000000 

0.000000000000 

0.000000000000 

o.oooooooooooo 

0.000000000000 

0.000000000000 

0.000000000000 


TABLE  IV 
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5.2.8  Sensitivity  and  Stability  Analysis 

For  a  small  positive  number  a  we  replace  the  problem  (5.19)  in  L\ s(0, 1)  with  homogeneous 
Neumann  boundary  conditions  by  the  following  problem  with  non-homogeneous  Neumann 
boundary  conditions: 

wt  -  twxx  +  wwx  =  0,  x  €  (0, 1),  (5.63) 

16(0,1),  t>  0 

wx(0,t)  =  wx(l,t)  =  -a,  a  >  0, 

w(x,  0)  =  </>(x). 

We  show  below  that  solutions  of  (5.63)  are  highly  sensitive  to  the  boundary  condition 

parameter  a.  „ 

For  4>  €  L2as(0,  1)  the  solution  w(-,t)  remains  in  L\s{ 0,1)  for  all  t  >  0  by  Theorem  5.3. 

Therefore,  the  Burgers’  problem  (5.63)  is  equivalent  to 

wt  -  ewxx  +  wwx  =  0,  x  G  (0, 1/2),  t  >  0  (5.64) 

wx(0,t)  =  -a,  w(l/2,t)  =  0,  a  >  0, 
w(x,  0)  =  <f>(x). 

The  stationary  problem  associated  with  (5.64)  isgiven  earlier  in  (5.30),  (5.31).  We  also 
note  that  in  L^s(0, 1)  we  can  replace  this  problem  by  the  equivalent  problem 

£VXX  vvx  —  0,  X  G  (0, 1),  (5-65) 

t>*(0)  =  -a,  «(l/2)  =  0.  (5-66) 


which,  for  numerical  work,  is  more  tenable. 

As  we  have  noted,  the  function  (5.27)  satisfies  stationary  Burgers  equation  and  the 
derivative  at  x  =  0, 1  satisfies  (5.31)  (and  (5.66))  provided 

h'(0)  =  h'(  1)  =  -a.  (5-67) 

This  amounts  to  finding  c0  that  satisfies  the  equation 

2>sech2  (^)  =  a.  (5-68) 

In  the  space  L\s{0, 1)  there  are  exactly  two  solutions  of  (5.68)  for  a  small  enough. 
Namely,  there  exist  «  0  and  >  0  giving 

/  , — ”  \ 

(5.69) 

(5.70) 


h<(x)  =  y^tanh  ^^-(1/2  -  x)J 
h>(x)  =  y^tanh  (^^.{1/2  -  x)^j 
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and  these  functions  satisfy  the  nonhomogeneous  boundary  conditions 

MO)  =  -<*:  Ml)  =  ~a- 

To  see  that  there  are  exactly  two  such  values  of  cq  for  small  a,  let  s  =  ^/cq/ (2\/2e)  so 
that  equation  (5.68)  becomes 

s2sech2(s)  = 

The  function  f(s)  =  s2sech2(s)  has  a  critical  value  at  s0  «  1.2.  This  allows  us  to  conclude 
that  the  maximum  value  of  /  is  Me  =  8esosech2(s0).  From  the  graph  of  /  in  Figure  17,  it 
is  clear  that  this  maximum  imposes  a  smallness  constraint  on  a.  Namely,  in  order  for  the 
conditions  in  (5.67)  to  be  satisfied,  we  need 

a  <  Mt. 

For  fixed  e  and  a  sufficiently  small,  we  see  that  there  are  two  solutions  given  by  (5.69)  and 
(5.70)  and  both  functions  satisfy  the  conditions  in  (5.67). 

The  solution  /*<(•)  is  very  nearly  the  zero  function,  whereas  the  solution  /i>(  )  is  not 
usually  small. 


FIG.  17.  Graph  of  /,  for  e  =  .01  and  a  =  .007 


FIG.  18.  Graph  of  for  e  =  .01  and  cq  =  7.7287e  -  05 


FIG.  19.  Graph  of  h>(-),  for  e  =  .01  and  Co  =  .0072 

In  each  of  Figures  18  and  19  there  is  actually  two  functions  plotted,  h  computed  numer¬ 
ically  and  also  from  the  formula  (5.27). 

A  complete  analysis  of  the  mathematical  validity  of  these  stationary  solutions  for  Burg¬ 
ers’  equation  involves  a  careful  analysis  of  the  long  time  behavior  of  solutions  to  the  dy¬ 
namical  system 


wt  -  ewxx  +  wwx  -  fa 
wx(0,t)  =  Wx(l,t)  =  0, 
u-(x,0)  =  4>{x ), 

fa  =  a(8Q-S1)GH-1(0,l) 

Here,  6a  denotes  the  (5-function  concentrated  at  x  =  a.  and  by  H~1( 0,1)  we  denote  the 
dual  of  Jfx(0, 1)  which  consists  of  all  distributions  from  tf_1(R)  whose  support  belongs  to 
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[0, 1].  For  small  a  and  small  initial  conditions  <f>,  one  knows  that  global  in  time  existence  of 
solutions  to  the  above  system  and  the  existence  of  a  compact  local  attractor.  Unfortunately, 
for  larger  initial  conditions  these  results  do  not  apply  for  fa  as  above. 


5.2.9  Linearization  about  Numerical  Stationary  Solution 

To  determine  the  stability  properties  of  these  stationary  solutions  we  follow  the  development 
for  a  similar  problem  given  by  Kreiss.  Namely,  after  a  reduction  to  the  interval  (0, 1/2),  we 
consider  the  spectral  analysis  of  the  linearization  of  the  Burgers’  system  (5.63)  in  L\s{0, 1/2) 
about  the  function  h(-)  in  (5.27).  We  will  show  that  for  a  fixed  e,  the  first  eigenvalue  of  this 
linear  problem  is  negative  if  Co  is  small.  But  for  larger  co  this  eigenvalue  becomes  positive 
and  then  decreases  monotonically  to  zero.  The  remaining  eigenvalues  of  the  linearized 
problem  are  all  negative.  Thus  for  co  small  the  stationary  solution  is  stable  and  for  the  larger 
value  of  co  the  corresponding  stationary  solution  is  unstable.  However,  the  first  eigenvalue, 
which  is  positive,  is  so  small  that  the  dynamics  still  can  converge  to  a  nonconstant  stationary 
value. 

To  this  end,  let 

w  —  h  +  Sz,  z  €  -^as^)  1)> 

substitute  this  into  (5.63)  and  collect  terms  of  order  one  in  5  to  obtain 

zt  —  ewxx  +  ( hz)x  =  0,  (5-71) 

*€(0,1),  t>  0 

Zx(0,t)  =  zx(l,t)  =  0,  a>0, 

z(x,0)  =  z0(x), 

where  generally  we  assume  that  zq  is  small  in  L\s( 0, 1). 

We  can  replace  this  problem  with  the  numerically  stable  problem  on  the  half  interval 
(0, 1/2)  given  by 

zt  -  ewxx  +  ( hz)x  =  0,  (5.72) 

x  G  (0, 1),  £  >  0 

z*(0,t)  =  *(l/2,i)  =  0,  a  >  0, 

z(x,  0)  =  z0(x). 


Associated  with  this  problem  we  consider  the  spectral  problem 

e<pxx  ~  {hf)x  =  -V,  <fix( 0)  =  ¥>(1/2)  =  0.  (5.73) 

It  is  useful  for  computations  to  replace  this  eigenvalue  problem  by  a  self-adjoint  problem, 
so  we  let 

f  rx  \  /  /oTT  / 1  \  \ 

(5.74) 


r)  =  exp 

and  seek  ip  in  the  form 


f  h(s) 
Jl/  2 


ds  =  sech 


ip  =  rjip. 
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After  some  straightforward  calculations  we  find  that  the  problem  (5.73)  can  be  replaced  by 
the  problem 

eipxx  -  #  =  AV>,  ^i(O)  =  ip(  1/2)  =  0.  (5.75) 


where 


Since  h  satisfies 


we  have 


9(x)  =  ^  +  ifc*. 

h? 

—hx  +  ^ —  co 


We  must  also  consider  the  transformation  of  the  boundary  conditions.  Note  that 


0  =  <px(0)  =  7,(0)  (0*(O)  +  — h(O)0(O)J  , 
0  =  v>(l/2)  =  77(l/2)^(l/2)  =  ^(1/2). 


Thus,  we  arrive  at  the  eigenvalue  problem 


Lfy'lp  =  0x,  ^7*0  A0  J 


subject  to  the  boundary  conditions 


where 


0,(0)  +  70(0)  =0,  0(1/2)  =  0 


7  =  ^-fc(O),  A  =  A/e,  p7  =  72(1  -  2sech2(7(l/2  -  *)). 

Zc 


As  Kreiss  noted  the  first  eigenvalue  of  (5.77)-(5.78)  satisfies 

.T.  ^  IIMJI 


(5.76) 


(5.77) 


(5.78) 


(5.79) 


for  all  functions  ip  G  L2( 0, 1/2)  that  are  sufficiently  smooth  and  satisfy  the  boundary  con- 
ditions  (5.78). 

Lemma  5.6  If  v/2co/(2e)  >  4,  then  there  are  •positive  constants  C  and  D,  independent  of 
e,  so  that  the  smallest  eigenvalue  of  (5.77)-(5. 78)  satisfies 


e 


(5.80) 


T](x)  -  T](x)  -  2t?(0) 

where  rj  is  defined  in  (5.74).  We  note  that  rj  is  a  symmetric  function  about  x  =  1/2,  it 
satisfies  the  boundary  conditions  (5.78)  and  is  smooth.  Also  we  have 

Lrj=  2qr)(0). 
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Thus,  we  have  the  estimate 


=  (j^\2v(0)q(x))2  dx\  ^  (j^\v(x)-2r1(0))2dx 


■( 


4cg 


< 


e2  cosh2(\/2co/(4e) 

C2 


(l  -  2sech2(\/2co/(2e)s))2  ds 


(  rl/2  _  2 

/  (2sech(\/2co/(4e))  -  sech(v^/(2e)s))  ds 

{Jo  J 


£2  exp  (2De) 


Here  we  have  used  the  fact  that  sech(\/2co/(2e)s))  is  a  monotone  decreasing  function  of 
s  and  the  assumption  7  =  \/2c^/(2e)  >  4  which  ensures  that 


(l-2sech2(7 s))2ds]  <  cq  (  [  '  (1  -  2sech2(7/2)sech(7s))2  ds 


Unlike  the  case  of  nonhomogeneous  Dirichlet  boundary  conditions  considered  in  by 
Kreiss,  for  nonhomogeneous  Neumann  conditions  the  resulting  linearization  about  the  equi¬ 
librium  is  not  exponentially  stable  for  all  7.  In  this  case  it  is  difficult  to  establish  this 
analytically,  so  we  only  present  the  numerical  results  obtained  for  the  first  three  eigenvalues 
computed  for  a  range  of  values  of  7  (he.,  <*>  for  fixed  e  =  .1).  For  one  value  of  0  we  have 
plotted  the  corresponding  normalized  eigenfunctions.  Note  that  when  cq  -  0  the  problem 
(5.77),  (5.78)  reduces  to 


</  =  \<p,  v'{  0)  =  0,  ¥>(1/2)  =  0 

with  eigenvalues  given  explicitly  by  \j  =  -(2 j  -  1)2tt2  for  j  =  1, 2,  •  •  •  and  with  associated 

eigenfunctions  q>j(x)  =  \/2cos((2 j  —  l)7rx)  for  j  =  1,2,  •  •  • . 

Thus,  for  small  7  the  problem  has  all  stable  eigenvalues.  However,  we  show  that  as  7 
increases  the  first  eigenvalue  3q,  as  a  function  of  7,  becomes  positive,  stays  positive  and 
decreases  to  zero  as  7  tends  to  infinity. 

We  have  fixed  e  =  .1,  set  7  =  and  varied  Co  from  co  =  0  to  Co  =  18  so  that  7 

varies  from  7  =  0  to  7  =  20.  The  eigenfunctions  are  depicted  in  Figure  20  with  7  =  16  (or 
cq  =  5.12)  while  Figure  21  contains  the  plot  of  7  versus  A j  for  j  =  1,2. 
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FIG.  20.  Eigenfunctions  with  with  7  =  16  (or  cq  —  5.12) 


FIG.  21.  Plot  of  Eigenvalues  7  versus  A j  for  j  =  1,2 


The  main  thing  we  learn  from  this  exercise  is  that  the  linearization  is  not  exponentially 
stable  for  large  values  of  Co  and  so  the  corresponding  stationary  solutions  are  not  stable  for 
the  Burgers’  system  with  nonhomogeneous  Neumann  boundary  conditions.  Nevertheless  for 
a  fixed  floating  point  accuracy,  as  we  have  seen,  the  numerical  solution  to  Burgers’  equation, 
for  certain  anti-symmetric  initial  data,  converge  to  the  solution  of  this  problem  for  some  cq. 


5*2,10  Conclusions 

In  this  paper  we  have  shown  that  problem  sensitivity  and  finite  precision  arithmetic  can 
combine  to  produce  false  numerical  solutions  to  steady  state  problems.  Although  one  might 
guess  that  it  is  possible  to  construct  pathological  examples  with  this  property,  it  is  somewhat 
remarkable  that  this  phenomenon  can  occur  for  a  “simple”  Burgers’  equation.  In  addition, 
in  the  space  of  antisymmetric  L 2  functions  L\ s(0, 1),  we  have  shown  that  the  steady  state 
Burgers’  equation  has  a  unique  solution  ( v(x )  =  0),  and4  yet,  discretized  versions  of  this 
equation  can  yield  non- unique  solutions.  More  importantly,  for  sufficiently  large  initial  data, 
marching  schemes  will  converge  to  this  discrete  (yet  false)  solution.  We  also  presented  a 
sensitivity  and  stability  analysis  that  provides  insight  into  these  numerical  difficulties. 

One  implication  of  these  results  is  that  more  analysis  needs  to  be  done  on  recent  “nu¬ 
merical  based  proofs”  of  nonuniqueness.  In  particular,  we  have  established  that  in  a  per¬ 
fectly  reasonable  model,  numerical  computations  that  yield  nonunique  discrete  stationary 
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solutions  do  not  imply  anything  about  nonuniqueness  for  the  continuous  boundary  value 
problem.  It  is  important  to  again  emphasize  that  it  is  finite  precision  arithmetic  that  causes 
these  false  solutions.  Therefore,  even  converging  an  algorithm  to  machine  roundoff  error 
will  not  eliminate  the  difficulty. 

Finally,  we  observe  that  the  false  numerical  solutions  can  differ  by  orders  of  magnitude 
from  real  stationary  solutions.  Therefore,  in  such  cases,  cascading  the  numerical  solutions 
into  an  optimization  or  control  algorithm  can  produce  bad  designs.  More  about  this  issue 
will  appear  in  future  papers. 
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M&A.E,  Cornell  University,  Ithaca,  NY  14853 


6.1  Background  and  Objectives 

The  Partnership  for  Research  Excellence  and  Transition  (PRET)  set  as  its  goal  to  produce 
outstanding  research  and  to  work  with  industrial  partners  and  Air  Force  partners  to  see  that 
the  research  gets  transitioned.  Because  of  the  nature  of  this  project,  the  personnel  working 
on  it  worked  closely  with  the  industrial  partners.  The  strategy  for  assuring  the  success 
of  the  transition  is  that  each  postdoc  was  assigned  as  a  “lead  postdoc”  for  one  industrial 
partner.  It’s  the  postdoc’s  job  to  assure  that  communication  and  research  transition  takes 
place. 

Our  proposal  stated  that  we  would  use  PDESolve  as  a  generic  tool  for  software  devel¬ 
opment.  BEAM  was  generous  to  provide  PDESolve  at  no  cost  for  use  by  PRET  members. 

The  goal  of  the  BEAM  -  Cornell  transitions  was  to  produce  PDESolve  modules  for  use 
by  the  PRET  team  and  for  applications  to  aerodynamic  design  and  materials  processing. 

6.2  Research  Achievements 

This  section  contains  a  summary  of  the  work  preformed  at  Cornell. 
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6.2.1  High  Level  3D  Navier-Stokes  Solver 

We  developed  a  3D  Navier  Stokes  solver  using  a  Galerkin  Least-Squares  to  circumvent 
the  usual  (LBB)  div-stability  restrictions  on  element  pairs.  A  finite  element  mesh  for 
this  problem  was  created  by  post-processing  a  model  generated  in  ProEngineer  (a  popular 
engineering  3D  CAD  package) .  A  number  of  computational  experiments  for  two  dimensional 
flows  were  compared  to  published  CFD  solutions  to  verify  the  results  of  this  software.  ,  and 
implemented  it  in  PDESolve. 

Because  of  high  level  symbolic  expression  specifications  of  PDESolve  the  entire  code 
was  less  than  200  lines  long.  This  is  very  impressive  because  a  similar  code  written  using 
traditional  techniques  such  as  FORTRAN  would  have  taken  three  years  or  more  to  write 
and  would  have  consisted  of  tens  of  thousands  of  lines  of  code. 

This  effort  pointed  out  a  deficiency  in  certain  aspects  of  PDESolve.  Specifically,  from 
this  effort  BEAM  learned  that  the  finite  element  engine  underlying  PDESolve  needs  to 
be  significantly  enhanced,  perhaps  even  re-written  if  it  is  to  support  state-of-the-art  DoD 
needs.  This  lesson  was  very  valuable  to  BEAM. 

This  effort  was  also  valuable  in  communicating  to  PRET  members  and  potential  DoD 
users  and  partners  the  power  of  mixed  symbolic-  numeric  computing  in  general  and  PDE¬ 
Solve  in  particular. 

As  part  of  the  Navier  Stokes  project,  an  evaluation  of  what  geometric  sensitivity  in¬ 
formation  would  need  to  be  provided  as  boundary  conditions  for  geometric  sensitivity  flow 
equations  was  done  as  well.  This  information  was  incorporated  in  the  design  of  a  2D  para¬ 
metric  geometry  engine  that  provides  sensitivity  information  as  well  as  interfaces  to  3D 
commercial  geometry  packages. 

6.2.2  2D  Shape  Optimization  and  Optimization  Architecture 

BEAM  Technologies  developed  a  parametric  2D  geometry  engine  that  provides  the  informa¬ 
tion  required  for  boundary  conditions  for  geometry  sensitivity  partial  differential  equations. 
We  applied  this  software  to  shape  optimization  model  problems.  The  results  were  excellent; 
we  were  able  to  implement  sensitivity  and  optimization  very  efficiently.  We  transitioned  the 
lessons  we  learned  to  BEAM  in  the  form  of  reccomendations  for  how  to  improve  PDESolve 
in  the  areas  of  optimization. 

Because  computing  cost  function  gradients  is  the  key  to  any  optimization  algorithm, 
we  developed  an  adaptive  finite  element  strategy  for  computing  cost  function  gradients  ef¬ 
ficiently  and  accurately.  Using  the  continuous  sensitivity  equation,  we  were  able  to  develop 
error  estimates  for  sensitivity  variables  (derivatives  of  flow  variables  with  respect  to  param¬ 
eters)  in  addition  to  error  estimates  for  the  flow.  Adaptive  mesh  refinement  was  performed 
to  reduce  both  sets  of  errors.  Thus,  we  were  able  to  compute  very  accurate  gradient  calcu¬ 
lations  using  fixed  computing  resources,  or  conversely,  computationally  efficient  algorithms 
to  obtain  a  desired  level  of  gradient  accuracy. 

At  the  time  of  the  research  project,  there  were  significant  discussions  in  the  community 
regarding  the  merits  of  the  analytic  sensitivities  compared  to  automatic  differentiation.  Au¬ 
tomatic  differentiation  has  the  advantage  of  being  applicable  to  existing  complex  scientific 
and  engineering  codes,  whereas  analytic  sensitivities  offer  improved  speed  (both  targeted 
at  computing  cost  function  gradients.)  To  bridge  this  divide  we  developed  a  method  for 
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Computation  of  continuous  sensitivity  equations  by  applying  standard  automatic  differenti¬ 
ation  software  to  (slightly)  modified  PDE  solvers  were  demonstrated.  Using  this  technique, 
we  were  able  to  use  AD  to  compute  derivatives  of  PDE  solutions  with  respect  to  shape 
parameters  without  the  need  to  compute  derivatives  of  discretization  parameters  (such  as 
the  mesh,  adaptive  time  stepping  algorithm,  etc.). 

6.2.3  Sensor  and  Actuator  Optimization  in  Flexible  Structure  Vibration  Con¬ 
trol 

An  area  of  importance  to  the  Air  Force  and  NASA  at  the  time  was  the  optimization  of 
the  location  of  sensors  and  actuators  for  smart  structure  control.  In  this  application  piezo¬ 
electric  sensors  and  actuators  are  placed  on  a  structure  and  a  control  algorithm  is  used  to 
dampen  the  vibrations.  This  allows  the  reduction  of  the  weight  of  a  sattelite,  for  example. 
Using  PDESolve  we  developed  innovative  technology  to  help  address  this  question.  This 
technology  was  transitioned  to  BEAM.  Research  included  algorithms  for  modeling  beams, 
plates,  box  beams,  and  three-dimensional  elastic  structures  with  distributed  SMEC’s  (Smart 
Materials  Embedded  Components)  for  controlling  vibration  and  stress.  Piezoelectric  actu¬ 
ators  are  modeled  as  three-dimensional,  homogeneous  bodies  using  the  quasi-static  linear 
theory  of  electroelasticity.  Governing  equations  for  the  substrate  structures  and  embedded 
components  are  expressed  in  finite-element  formulations.  The  software  employs  model  or¬ 
der  reduction  and  has  direct  interfaces  with  control  design  software  like  Matlab  and  CAD 
software  like  ProEngineer.  See  attached  report. 


6.2.4  Materials  Processing  Applications  of  Relevance  to  the  Air  Force 

Since  one  of  the  objectives  of  the  Cornell  -  BEAM  collaboration  was  to  develop  modules 
in  the  area  of  materials  processing,  an  effort  was  put  into  identifying  materials  processing 
applications  that  may  be  of  interest  to  BEAM  on  the  one  hand  and  are  relevant  to  the  Air 
Force  on  the  other  hand.  Two  applications  were  identified: 

Injection  molding  of  viscoelastic  and  viscoplastic  materials  As  a  result  of  collabora¬ 
tion  with  Allied  Signal  we  learned  that  Allied  Signal  has  developed  a  certain  ceramic  slurry 
that  they  were  interested  in  injection-molding  into  shapes  that  would  eventually  become 
turbine  blades  after  a  baking  process.  Ceramic  blades  was  an  area  DARPA  made  heavy  in¬ 
vestments  in,  and  Allied  Signal  was  one  of  the  beneficiaries  of  these  funds.  However,  despite 
the  success  in  the  development  of  the  basic  materials,  Allied  Signal  had  a  problem  with  the 
manufacturing  process.  Specifically,  optimization  of  the  location  of  the  material  injection 
point  was  an  issue.  A  collaborative  effort  was  put  in  place  to  try  and  solve  a  generic  3D 
injection  molding  problem,  followed  by  models  with  more  complex  material  physics.  This 
effort  did  not  yield  conclusive  results,  as  PDESolve  had  difficulties  in  converging  some  o 

the  problems.  _ 

Chemical  Vapor  Infiltration  in  the  Manufacturing  of  Carbon-Carbon  Compos- 

ites 

Carbon-Carbon  composites  is  a  technology  the  Air  Force  has  made  significant  invest¬ 
ments  in.  One  of  the  problems  with  this  technology  is  that  setting  manufacturing  parameters 
is  a  trial  and  error  process.  The  process  is  very  expensive  because  each  manufacturing  run 
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may  take  up  to  a  week  of  “baking  .”  The  objective  of  this  work  was  to  develop  a  capabil¬ 
ity  to  optimize  the  process.  The  first  step  was  to  set  up  a  simulation  of  the  process.  This 
project  was  a  great  success  -  the  symbolic  nature  of  PDESolve  enabled  solving  fully  implicit 
chemistry  using  a  species-sensitivity  matrix.  All  the  simulations  were  in  2D.  This  project 
again  proved  the  power  of  a  mixed  symbolic-numeric  approach.  This  was  a  collaborative 
effort  with  BEAM. 

6.2.5  Fluid  Structure  Interaction 

A  new  method  for  predicting  dynamic  effects  of  fluid-structure  interaction  was  developed 
by  BEAM,  for  the  purpose  of  optimizing  designs  to  minimize  flutter,  and  enable  the  control 
of  flutter  with  sensors  and  actuators.  BEAM  made  several  extensions  to  PDESolve,  as  well 
as  a  PDESolve  code  generator  using  Mathematica.  These  were  provided  to  Cornell  at  no 
cost  for  the  purpose  of  researching  BEAM’S  new  method  and  providing  feedback  to  BEAM. 
The  new  method  holds  great  promise,  and  is  currently  the  only  known  approach  that  can 
account  for  viscous  effects  on  flutter  and  limit  cycle  oscillations  using  CFD-quality  fluid 
dynamics  data.  The  results  of  this  research  are  summarized  in  the  papers  by  Miller. 

6.2.6  Control  and  fluid  Structure  interaction 

In  work  related  to  both  flexible  structure  vibration  control  and  to  fluid  structure  interaction 
(see  sections  6.2.3,  6.2.5)  we  examined  the  structure  of  the  turbulent  boundary  layer  with 
an  eye  to  its  control  and  modification  over  flexible  structures.  As  a  first  step,  we  considered 
a  compliant  surface.  Specifically,  we  followed  an  approach  in  which  only  the  large  scales 
of  the  boundary  layer  are  resolved,  the  smaller  scales  being  parameterized.  In  this  way,  a 
low-dimensional  model  of  the  layer  is  developed,  which  can  be  investigated  dynamically  as 
it  interacts  with  various  structures.  The  general  approach  is  described  in  the  papers  below. 

More  general  contributions  are  contained  in  the  papers  by  Lumley  and  his  co-workers. 
In  particular,  the  book  by  Lumley  describes  the  application  of  some  of  these  same  ideas 
to  flow  in  engine  cylinders.  A  tumbling  flow  leads  to  elliptic  instability  which  promotes 
high  turbulence,  increasing  flame  speed  and  reducing  pollution.  Our  approach  permits  the 
examination  of  this  dynamical  process  and  its  manipulation  to  optimize  these  effects. 

The  second  book  by  Lumley  is  a  report  on  current  interesting  research  directions  :  in  fluid 
dynamics  prepared  under  the  auspices  of  the  U.  S.  National  Committee  on  Theoretical  and 
Applied  Mechanics.  This  was  intended  for  members  of  Congress  and  their  staffs. 
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7  Honors 

•  Professor  John  Burns  was  named  Hatcher  Professor, 

•  Professor  John  Burns  was  elected  Fellow  of  the  IEEE, 

•  Professor  Jeff  Borggaard  was  named  a  PECASE  Fellow, 

8  Personnel  Supported 

The  following  people  were  supported  in  part  under  Grant  F49620-96-1-0329: 
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9  Interactions/Transitions 

Our  efforts  to  expedite  the  transition  of  our  research  to  industrial  and  Air  Force  needs  are 
manifested  by  direct  industrial/laboratory  interactions  and  participation  at  professional 
meetings.  One  of  the  major  components  of  this  effort  was  active  cooperation  and  coordina¬ 
tion  with  the  Air  Force  Research  Laboratory  (AFRL)  and  with  our  industrial  partners.  We 
have  actively  worked  with  all  of  our  industrial  partners  and  with  groups  at  AFRL/VACA 
and  AFRL/DE. 
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Industry /Laboratory  Interactions 

Aerosoft  Inc, 

We  worked  closely  with  AeroSoft  on  several  projects  to  develop  software  packages  for  anal¬ 
ysis  and  design  of  aerospace  systems.  The  central  theme  in  these  efforts  is  the  continu¬ 
ous  sensitivity  equation  method  (SEM)  for  approximating  the  effect  of  parametric  design 
changes  on  aerodynamic  performance.  As  of  this  writing  the  production  code  SENSE  is 
at  Version  1.0.4.  Recent  additions  permit  sensitivity  calculations  in  turbulent  flows  with 
one  and  two-equation  turbulence  models.  We  continue  to  support  AeroSoft’s  work  with 
AFRL/DE  on  design  and  analysis  for  COIL  lasers.  As  noted  above,  we  are  studying  formu¬ 
lations  to  efficiently  couple  single  discipline  sensitivity  codes  for  the  study  of  multi-physics 
problems.  This  work  is  being  supported  by  AFRL/DE  and  is  described  below. 

Boeing  Defense  and  Space  Group  (BDSG) 

Recent  efforts  with  our  Boeing  partners  are  focused  on  transition  of  our  recent  work  on 
CFD/Sensitivity  methods  for  estimating  rotary  aerodynamic  derivatives.  Dr.  A.C.  Li- 
mache  has  completed  his  Ph.D.  studies  in  this  area  and  research  is  being  transitioned  to 
Boeing.  Our  initial  objective  is  to  implement  required  changes  in  a  3D  Euler  code  to  pro¬ 
vide  a  capability  to  estimate  aerodynamic  forces  and  moments  in  a  generalized  steady  flight 
maneuver. 

Directed  Energy  Directorate  of  AFRL 

Dr.  T.  J.  Madden  and  others  at  (AFRL/DELC)  are  involved  in  efforts  to  develop  technolo¬ 
gies  for  improved  performance  in  chemical  oxygen-iodine  lasers  (COIL).  Gaseous  chemical 
lasers  can  provide  lightweight,  efficient  energy  sources  for  a  wide  variety  of  Air  Force  systems 
including  airborne  and  space-based  directed  energy  weapons.  Computer-based  design  tools 
can  lead  to  rapid  development  of  efficient  laser-power  systems.  In  addition  to  the  coupled 
sensitivity  analyses  noted  above,  we  are  developing  an  alternative  formulation,  based  on  a 
paraxial  wave  equation,  for  modeling  energy  extraction  in  the  laser  cavity.  This  replaces 
a  “discrete”  ray-tracing  algorithm  currently  in  use  at  AFRL/DELC  and  is  a  more  natural 
setting  for  continuous  sensitivity  analysis. 

Air  Vehicles  Directorate  of  AFRL 

Dr.  Siva  Banda  and  others  at  AFRL /VAC A  are  starting  a  new  effort  to  develop  control 
technologies  for  application  to  flow  control.  We  are  working  with  Dr.  Banda’s  group  on 
flow  control  and  computational  tools  for  design  of  distributed  parameter  systems.  We  plan 
to  extend  earlier  work  on  functional  gain  computations  to  a  practical  experimental  test. 
Dr.  Burns  will  be  spending  time  at  AFRL/VACA  to  help  initiate  this  project. 

10  Publications  Produced  Under  the  Grant 


Books 

1.  Computational  Methods  for  Optimal  Design  and  Control,  Edited  by  Jeff  Borggaard, 
John  Burns.  Eugene  Cliff  and  Scott  Schreck,  Progress  in  Systems  and  Control  Theory, 
Birkhauser.  Boston,  1998,  475  pages. 
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